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SUMMARY 


Compared to numerical analysis, symbolic and algebraic manipulation is unfamiliar to 
people in the research field. As its name implies, symbolic and algebraic manipulation can be 
simply interpreted as a computerized operation which can retain symbols throughout 
computations and express results in terms of symbolic forms. For example, the coefficients a, 
b, and c in the quadratic polynomial equation ax*+bx+c=Odo not need to be known in order to 
find its roots. The equation itself can be directly input to a computer and the results will be 

(-b + y/b‘-4ac) . ( - b - V b‘ - 4ac 1 

x < 3 and 23 

If the numerical values are required, three coefficients can be specified and solutions will be 
expressed as numbers. 

From the example above, at least two unique characteristics of symbolic and algebraic 
manipulation can be observed. First unlike numerical analysis, the solutions from symbolic and 
algebraic manipulation are exact and therefore no round-off error is introduced. Second, the 
solutions are the same as those derived by hand. Therefore the extension of human capability to 
handle more sophisticated fomulations becomes feasible by computer. 

In the first chapter of this report, the history of symbolic and algebraic manipulation is 
introduced. The sencond chapter chronologically reviews the literature regarding the application 
of symbolic and algebraic manipulation in the engineering field. The capabilities of symbolic 
and algebraic manipulators are demonstrated in chapter three by selected examples. Chapters 
four through six demonstrate applications of symbolic and algebraic manipulation. Chapter 
four describes the automatic formulation of applied mechanics problems, chapter five covers 
the materially nonlinear, rigid-plastic ring compression problem, and chapter six discusses 
plate problems. The final chapter summarizes the overall conclusions of this report. 

It is well known that there are some difficuilties existing in the symbolic and algebraic 
field. The report proposes a remedy to avoid the difficulties and successfully accomplishes the 
applications. Due to this breakthrough, the solution of some previously insol vable problems 
become available. In addition, one of the advantages found in this research is believed to be 
crucial for improving the execution efficiency of numerical programs. 



CHAPTER I 


HISTORY OF SYMBOLIC AND ALGEBRAIC MANIPULATION 


1.1 Introduction 

Symbolic and Algebraic Manipulation (abbreviated as SAM) software is one of the new 
products of modem technology for use with highly developed digital computers. Traditionalists 
might say that SAM software is a misuse of modem computers. This is true if the viewpoint is 
adopted that a computer is a machine which only counts numbers. This viewpoint, however, 
severely limits the emerging artificial intelligence capabilities of computers. For instance, in 
addition to numbers, there are many symbols which define appropriate mathematical relations 
in a calculus book. Can we ask a computer to do these analytical derivations for us? This is a 
great question which finally led to the birth of SAM and added "soul" to the computer, to make 
it think more like a human brain. This idea, which originated before 1953, has had an impact 
on a variety of fields, such as science, industry and education. Therefore, although its history 

is not as long as that of the classical sciences, its impact has been so large that a record of its 
history is deserving. 

At the initiation of this dissertation effort in 1986, there were already a number of SAM 
systems available on the market. Some of them were more than ten years old, such as 
FORMAC, REDUCE and MACSYMA. Others were just being developed, such as muMATH 
and MATHEMATICA. Ironically, most relevant documents either ignore the history of these 
systems or just skim over them briefly. Only one book, written in 1969 by Jean E. Sammet 
[1], includes historical details, however, it is too old to cover recent developments. Most of the 
major systems used today have been produced since then. Therefore, it is necessary to collect, 
rewrite, and update the history of the SAM systems. It is hoped that the interested researchers 
will get a complete picture of the development of SAM systems. Through an understanding of 
the history they will be able to grasp the direction of the field and devote themselves towards 
making a further contribution. This is the major purpose of the chapter. It is much more 
important than just knowing how to run the SAM systems. 
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The first idea for using computers to do SAM can be traced back to two master theses 
published in 1953 ([2], [3]). Three years later, what is believed to be the earliest SAM system 
called PM, was developed at IBM [4], Now there are many SAM systems on the market for 
various computers. Some of them are designed for general purpose usage, while others have 
been developed for particular applications. Generally speaking, the evolution of symbolic and 
algebraic manipulation can be classified into three stages. They are : 

1. The first generation (1953-1965)— software to appear in this generation was PM, 

ALGY, FORMAC, MATHLAB and ALTRAN. Because of the limitation of hardware 
capacity, the systems in this stage were small in size and immature in content. 
Therefore, most of them became obsolete or were revised. 

2. The 2nd generation (1966-1975)— -software to appear in this generation was REDUCE 

and MACSYMA. These systems took advantage of the improvement of hardware 
memory capacity. They contain many built-in functions, are large and are for general 
purpose usage. All ol them run on the mainframe. 

3. The 3rd generation (1976- present)— -some representatives of more recent systems are 

muMATH, MATHEMATICA, and DERIVE. Unlike the systems of the second 
generation, the systems in this generation are designed to run on microcomputers. This 
has been possible due not only to the improvement of memory capacity in 
microcomputers, but also due to the requirement of most users who just need quick 
checks or moderate manipulations. 


Details of the histories of the systems will be described in the following subsections 
individually. For the sake of clarification, a summary is also included in Table 1.1. 

1.2 History of SAM systems 

1.2.1 PM 

PM is believed to be the earliest computerized algebraic system in the world. It was 
developed by George E. Collins at the IBM research center in Yorktown Heights, New York. 
Although the first document was published in 1966 [4] its beginning dates back to 1956. 
Written in assembly language for the IBM 701 computer, PM contained the subroutines for 
addition, subtraction and multiplication of multiple-precision integers, and subroutines for 
performing the same operations on multivariate polynomials with multiple-precision integer 
coefficients. Between 1956 and 1966 PM was reprogrammed for the newer IBM computers 
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(e.g. 709, 7090 and 7094), and augmented to include new operations (such as integer gcd) 
with various improvements (e.g. the incorporations of list processing and dynamic storage 
allocation). In 1966 Dr. Collins became a professor of computer science at the University of 
Wisconsin. With the aid of graduate students, the PM system was converted to the SAC-1 
system in 1973. In spite of its eventual replacement, PM, as the first SAM system, was still 
very significant in the field. 

1.2.2 ALGY 

Although ALGY has few functions, it was one of the earliest SAM systems in the 
world. The developmental work was started at Western Development Lab-Philco Co. in Palo 
Alto, California by Bemick, Callender and Sanford, around 1961 [5], It was interactive and 
allowed expressions written in a notation similar to FORTRAN as input, with some deviations. 
For example, the $ was used instead of ** to represent exponentiation, and all natural integers 
were expressed as fractions, for example, 0 and 1 were denoted by 0/1 and 1/1, respectively. It 
only contained a few commands, such as : 

• OPEN : expanding the expression in the parenthesis 

• SBST : making substitution 

• FCTR : factoring a given expression 

• TRGA : expanding th esin(a+b) into sin(a)cos(b)+sin(b)cos(a) 

Which is why the authors said that only two hours instruction was enough to use it. 
Although ALGY didn't come into extensive use, some of its ideas were succeeded by the 
FORMAC system, which is still popular today. 

1.2.3 FORMAC 

FORMAC is an acronym of FORmula MAnipulation Compiler. It was developed by J. 
E. Sammet and Robert G. Tobey at IBM's Boston Advanced Programming Department in 
July, 1962 [1], Five months later (December, 1962), the first complete draft of language 
specifications was prepared and implementation design started immediately thereafter. After 18 
months of extensive experiments, the first complete version was successfully running on the 
IBM 7090/94 computer in April, 1964. For the sake of obtaining feedback from users to make 
further improvement to the system, FORMAC was released for public use by the authors 
themselves (not by IBM) in November, 1964. This version of FORMAC was written in 
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FORTRAN IV. Three years later (November, 1967), the new version of FORMAC written in 
PL/I was released by the authors for use on IBM/360 systems. 

The FORTRAN version of FORMAC kept most commands and notations of 
FORTRAN IV. In addition, there were a couple of new commands added to allow it to do 
algebraic manipulations. For example, LET assigns symbols to variables instead of numbers in 
FORTRAN, SUBST makes substitution, EXPAND removes all parentheses in expressions 
and COEFF obtains the coefficients of variables. The major PL/I FORMAC capabilities can be 
divided into the following categories : 

1. User control of simplification : EXPAND for expanding the parentheses expression, 
DIST for applying the distributive law to all products of sums, etc. 

2. Substitution : EVAL(expr,a,b) replaces a in expr by b. 

3. Differentiation : DERIV performs partial differentiation. 

4. Expression analysis : COEFF(exprl,expr2) returns the coefficient of exprl in expr2. 
NUM and DENOM return the numerator and denominator, respectively. HIGHPOW 
and LOWPOW return the highest and lowest power. 

5. Storage allocation : SAVE(var) for storing the var to secondary storage. 

6. Output : PRINT_OUT(expr) to print out the required expressions. 

7. Built-in functions : these include trigonometric, logarithm, exponentiation, square root, 
hyperbolic function, etc. 

8. User defined function : the user can define functions as needed. 

FORMAC is now one of the most popular SAM systems. It is the first reasonably 
general purpose system to receive extensive usage worldwide. With the advantages of longer 
history and larger numbers of users, its accumulated contributions to SAM field are 
remarkable. In 1977, the new version, called FORMAC 73, was released to replace the old 
one. 
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1.2.4 MATHLAB 


MATHLAB 1 system was developed by C. Engelman and his employees at MITRE Co 
m 1964 [6], Its source language is LISP, but the commands am defined as English words. For 

ls the co "’ raa ' ,d ,0 simpl,fy x and name T- in the fall 
’ e lrst version °f MATHLAB was replaced by the second version, MATHLAB 68 

whtch operated on a PDP-6 machine with 256 K core memory. The input and output were 

h rough a teletype-hire keyboard wtlh a fixed character display scope. The notations m the 

second vers, on were more ALGOL-like. MATHLAB was the fits, complete on-line system. 

1.2.5 ALTRAN 

ALTRAN is a system developed at the BELL TELEPHONE Laboratory in Murray Hill 
New Jersey by W. S. Brown, M. D. Mcllroy, D. C. Leagues and G. S. Stoller [1] „ was 

™”‘"® H 964 °" lhe ,BM 7090/7094 ' 7040/44 ’ “ The baSIC '“"Suages whtch 
ALTRAN adapted were a mixture of FORTRAN II and FORTRAN IV. Smce it was hmited to 

use in the BELL Lab., its contributions to the SAM field were small. 

1.2.6 REDUCE 

REDUCE was developed by A. C. Hearn of Rand Corporation, Caltfomia, ,n 1963 

,. “™’ he me ‘ Dr - John McCanh >'- “ of Uk LISP language, who suggested 

e use of LISP for the problems of elementary particle physics. Since then. Dr. Hearn as a 
theoretical phystest, has worked ,n the SAM area. In August 1966, the firs, publication was 
issued [81. This paper only talked about the spectfic appl, cation of SAM techniques to 
elementary parade phys.es. Two years later (1968), the first paper describing a general algebra 
system, -REDUCE", was published [9J. The name of REDUCE onginated from thts paper Its 
name ,s no, an acronym. According to the description from the author himself, its name was 
actually ended as a wit. He satd -algebra system then as now,tended to produce very large 
expressions for many problems, rather then reduce the results to a more manageable form- 
The system at thts „me was called REDUCE for distinction from the new vetsion REDUCE •> 
w tch apjKamd ,970. The big improvement was tha, the whole system was written a'n 
t e talect (call RLISP), rather than the parenthesized notation of LISP in which 
REDUCE was written. A, thts time, the REDUCE 2 system was also released to users, making 
the beginnings of a user community. Thereafter, REDUCE 2 was implemented successfully on 
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the Michigan Terminal System (MTS) of the University of Michigan by Mike Alexander. After 
a long silence, REDUCE 3 was distributed in 1983. Several significantly new packages were 
added in this version, such as analytic integration, multivariate factorization, arbitrary precision 
real arithmetic and equation solving. Following REDUCE 3, upgraded versions were also 
released. They were REDUCE 3.1 released in 1984, REDUCE 3.2 in April, 1985, REDUCE 
3.3 on July 15, 1987. Each of them contains bug fixes and additional capabilities. Instead of 
implementation on MTS, the REDUCE 3.3 was first implemented on the APOLLO workstation 
in the Computer Aided Engineering Network (CAEN) of the University of Michigan. 
REDUCE 3.3 was also updated once in January 15, 1988. 

REDUCE system has become one of the most well-known SAM systems. Its general 
purpose design makes it possible to be used in a wide variety of areas. Its contributions are 
confirmed by the number of papers published in different fields. 

1.2.7 SCHOONSCHIP 

SCHOONSCHIP was designed by M. Veltman at CERN, Switzerland in 1964 [10]. 
Its major applications are in the field of high energy physics, but it is sufficiently general to be 
used for other calculations. It can deal easily with expressions of 10 4 to 10 5 terms on the CDC 
6000 computer. It was limited to use within CERN. 

1.2.8 ANALITIK 

ANALITIK was developed at the Institute of Cybernetics in Kiev, Soviet Union, by the 
direction of the well known Soviet cybernetician and academician V. M. Glushkov [11], The 
first paper discussing the system features was published in 1964. The language it used was 
ALGOL-like and close to that of traditional mathematical notation and natural language. It 
possessed interactive and batch processing modes. Since its implementation is highly machine 
dependent, ANALITIK has only run on the MIR-2 computer. 

1.2.9 FLAP 

FLAP was written in LISP 1.5 by A. H. Morris, Jr. at the U.S. Naval Weapons 
Laboratory in Dahlgren, VA. prior to 1967 [1]. Obviously, the FLAP system wasn't released 
to the public. 
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1.2.10 SAC 


The SAC system was developed by Dr. George E. Collins at the University of 
Wisconsin, Madison. The first version, SAC-1, was distributed in 1967 [4], This was a highly 
portable general purpose system, developed to replace one of the very earliest computer algebra 
system, PM in IBM [see 1.2.1], SAC-1 was replaced by SAC-2 in July, 1980. The SAC-2 
was programmed in ALDES language, which was designed by Rudiger Loos and G.E. Collins 
in 1973 to 1974. The SAC-2 system also provided the translator from ALDES to standard 
FORTRAN to maintain its portability. 

1.2.11 MACSYMA 

MACSYMA is an acronym of project MAC's SYmbolic MAmpulator. It was originally 
designed by C. Engelman, W. Matin, J. Moses for project MAC at M.I.T. in 1968. The 
implementation of it began in July, 1969. The system has quintupled in size since the first 
paper describing it appeared in 1971 ([12] [13] [14] [15]). It was made available over the 
ARPA networks in May, 1972. MACSYMA has a lot of built-in mathematical functions and 
graphic facilities which have made it one of the most powerful SAM systems in the world. 
Unfortunately, the University of Michigan didn't have it until September, 1988. The one 
implemented on the APOLLO workstation in CAEN of the University of Michigan still doesn't 
have a graphics package. 

1.2.12 SCRATCHPAD 

Although the name of SCRATCHPAD was chosen in 1970, the initial work on it can 
be traced back to 1965. The SCRATCHPAD system was designed principally by James H. 
Griesmer, Richard D. Jenks, Fred Blair, David Yun, and their colleagues, at the IBM Thomas 
J. Watson Research Center, Yorktown Heights, New York ([16] [17] [18]). Unfortunately, 
the name of SCRATCHPAD was not used for the first paper, presented in Bonn in 1970. One 
year later, a revised version by Dick Jenks, called SCRATHPAD/1, was demonstrated at 
SYMSAM/II in March 1971. After combining some new features, such as history file 
(allowing users to backtrack), and system commands, the first completed SCRATCHPAD/ 1 
manual was eventually published in 1975. After this, there seemed a stagnation in the progress 
of the SCRATCHPAD system due to personnel changes. Jim Griesmer left the group to be a 
manager of education at IBM research and Dick Jenks went to the University of Utah for a 
sabbatical. When Dick Jenks returned to Yorktown Heights in the fall of 1977, David Yun 
agreed to organize the "mode-base" ideas originated by Dick Jenks in 1973. This led to the 
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NEWSPAD which thereafter was renamed to SCRATCHPAD 84 at the New York conference 
in 1984. However, the name of SCRATCHPAD 84 was not quite appropriate since it would 
take more than one year to finish the system. Therefore it was changed into SCRATCHPAD II, 
which is the name used now. It became available in 1985 for test and evaluation to a limited 
number of users from an IBM owned mainframe via telenet, CSNET and ARPANET. As yet, 
it is not commercially available. 

1.2.13 CAMAC 

The CAMAC system was designed by Vera Pless in 1973 at M.I.T. [19]. The first 
version of it ran interactively and was written in FORTRAN with sections in assembler 
language. When Vera Pless moved to Chicago Circle in 1975, the CAMAC system was 
transferred to the Circle's IBM 370-158 by William Pattern. The name of CAMAC is an 
acronym of Combinatorial and Algebraic Machine Aided Computation. As the name implies, it 
was for a specific application. 

1.2.14 SHEEP 

SHEEP was designed by I. Fnck at the University of Stockholm, Sweden in 1975 
[(20], [21]) It was specialized for manipulating components of tensors. The source language it 
uses is MACRO- 10. It runs on DEC PDP 10 and PDP 20. The first version of SHEEP, now 
called SHEEP 1, was written in assembler code for the DEC- 10/20 computer. Unlike the first 
version, SHEEP 2 is written in standard LISP. 

1.2.15 ORTOCARTAN 

ORTOCARTAN is written in LISP. It was designed by Andrzej Krasinski in Poland in 
1977 [22], Its name is an acronym and is due to the specific application to the calculation of 
Riemann, Ricci, Einstein and Weyl tensors from a given metric tensor using an ORTHonormal 
set of CARTAN forms. Although the author said it could be relatively easily extended for other 
uses, such as inverting matrices of arbitrary rank, it did not come into wide use. 

1.2.16 MAPLE 

The MAPLE system was designed by Bruce Char, Keith Geddes, W. Morven 
Gentleman and Gaston Gonnet at University of Waterloo, Canada in December 1980 ([23] 
[24]). The name "MAPLE" is not an acronym but rather it was simply chosen as a name with a 
Canadian identity. There were two goals which oriented MAPLE's design. The first was to be 
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used on a time sharing mainframe computer. The second was to run it on a microprocessor- 
based workstation. This was the major difference between MAPLE and REDUCE (or 
MACSYMA). 

1.2.17 muMATH 

Written in muSIMP (a LISP-like language), muMATH was designed and developed by 
David R. Stoutmeyer and Albert Rich at the University of Hawaii in 1977 [25], The first 
version was called muMATH-77 and was experimental. Two years later, the Software House, 
Inc., was founded by David Stoutmeyer and the first product, muMATH-79, was distributed 
to users for the CP/M-80 operation system or Apple II family Z80 machine with 64 K bytes 
core memory required. The second product called, muMATH-83, was not released for the IBM 
personal computer until 1983. The muMATH-83 needs 256 K bytes RAM memory. Recently, 
DERIVE has taken over the place of muMATH-83. The significant improvement is that 
DERIVE combines the numerical, algebraic and graphical functions together, rather than just 
algebraic functions of muMATH. The DERIVE system requires 512 K memory space for 
normal execution. 

1.2.18 MATHLIB & SMP 

MATHLIB is an interactive general purpose SAM system. It was originally designed 
and developed under the auspices of the Department of Mathematics at Harvey Mudd College, 
California, in 1978. It was one of the products of Innosoft International Inc., of Claremont, 
California and became commercially available in 1983. It can perform numerical and symbolic 
operations. In addition, its graphical output is device-independent and allows it to be processed 
by over 150 different graphics devices. The PRS subroutine embodied in the algebraic 
subsystem of MATHLIB has more than 250 built-in functions for manipulation of 
mathematical expressions. 

SMP is another product of Innosoft International Inc.. It was designed at the California 
Institute of Technology. It's written in C language and was originally developed to run on 
VAX/780 under the UNIX operating system. In addition, there is a special design character in 
SMP to allow for easy conversion between operating systems. It needs at least 2.5 megabytes 
of memory space for typical usage. 
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1.2.19 MATHEMATICA 


MATHEMATICA is a recent product of Wolfram Research, Inc. There are several 
versions of MATHEMATICA for a variety of computers, such as for Apple Macintosh, DEC 
VAX, IBM ,Cray, and so forth. It was designed and implemented by Stephen Wolfram, Daniel 
Grayson, Roman E. Maeder, and their colleagues at the University of Illinois in 1988 [26], It 
integrates the algebraic manipulation, numerical computation, and graphical functions together 
and allows the resultant expressions to be outputed in C code, FORTRAN code, and text form. 
Its source language is C. The memory requirement for normal operation is about 3.7 mega 
bytes. The MATHEMATICA as well as DERIVE are expected to be two dominant systems in 
the coming decade. 

1.2.20 Mathcad 

Mathcad is developed by Mathsoft, Inc. at Cambridge, Massachusetts. The earlier 
version appeared on market around 1987. This system adopted the core functions of MAPLE 
and extented itself by including the graphic capacity. It is written in C language. The latest 
version 3. 1 is available in 1992. This newest version can run in IBM, Macintosh PCs and unix 
based machines. The minimum space requirements are two megabyte RAM and seven 
megabyte hard disk. Unlike most of the SAM systems, the command inputs in this system are 
menu driven. This allows users to communicate with machine by simply picking and clicking. 
This unique feature not only saves users lots of efforts in typing but also reduces human errors 
which sometimes turn out a unmanageable, hard- to-be-debugged results. 

1.3 Conclusion 

From a history of the SAM system, we can draw the following conclusions : 

• The SAM systems evolved from small, immature systems to well-designed, multi- 
function systems, to compact systems which can be used on microcomputers. 

• At present, there is no unique best system. The definition of the best SAM system 
depends on many variables, such as computer availability, availability of software, 
familiarization with software, the problem to be solved, software contents, circumference 
facility, and so forth. 
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system 


year 


remarks 


PM 

1956 

IBM 

ALGY 

1961 

WDLP Co. 

FORMAC 

1962 

IBM - Boston 

MATHLAB 

1964 

MITRE Co. 

ALTRAN 

1964 

BELL Lab. 

REDUCE 

1963 

Rand Co. 

SCHOONCHIP 

1964 

CERN 

ANALITIK 

1964 

Soviet Unions 

FLAP 

1967 

U.S. Navy 

SAC 

1967 

Uni. of Wisconsin 

MACSYMA 

1968 

M.I.T. 

SCRATCHPAD 

1965 

IBM-Yorktown Heights 

CAMAC 

1973 

Vera Pless 

SHEEP 

1975 

Sweden 

ORTOCARTAN 

1977 

Poland 

MAPLE 

1980 

Canada 

muMATH & DERIVE 

1977 

Uni. of Hawaii 

MATHLIB 

1977 

Harvey Mudd College 

SMP 

1977 

Caltech 

MATHEMATICA 

1988 

Uni. of Illinois 

Mathcad 3. 1 

1992 

Mathsoft Inc. 


Table 1.1: List of symbolic and algebraic systems 
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CHAPTER II 


SURVEY OF THE LITERATURE ON SYMBOLIC AND ALGEBRAIC 

MANIPULATION 


2.1 Introduction 

The documents published in the symbolic and algebraic manipulation field are not as 
plentiful as those in the area of numerical analysis. However, after a careful classification of the 
existing documents, one finds that the developmental history is closely related to research 
directions and content of the publications. In general, the documents about SAM may be 
divided into four categories. They are : 

1. About SAM system itself — More than half of the existing papers belong to this class. 
Most of them were published in the period of the first generation. The contents are 
focused on the following topics : 

(a) The introduction of the new SAM system, including the capacities, functions, etc. 

[1][2][3], 

(b) The technical reports of softwares [4] [5] [6] [7] . 

(c) The data structure, language and implementation [8][9], 

2. Applications to science — This class of publications is the second largest of the existing 
SAM papers. One of the major impetuses in developing SAM systems was due to the 
requirements from scientists, especially in the fields of elementary particle, general 
relativity and celestial mechanics. Some of the famous examples were collected in the 
paper by Hearn [10]. One of them is the recalculation of Delaunay's moon coordinates 
by Deprit, Henrard and Rom in 1970 [11] The others are such as Campbell and Hearn's 
analysis of the Feyman diagram [13], Rudiger Loos' work about Archimedes' cattle 
problem [14], and Roberts' and Boris' on the solution of partial differential equations 
[15]. 
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3. Applications to engineering — In fact, most engineering problems are not solvable 
analytically. Therefore, numerical approximation usually predominates in the solution of 
engineering problems. This is one of the reasons why the engineering applications of 
SAM has not been as popular as those on science. However, there are two factors which 
necessitate the use of symbolic and algebraic manipulation in engineering. The first is 
that the accuracy requirement of a solution of numerical approximation becomes more 
and more strict today. The second is that the problem to be solved usually involves more 
sophisticated algebraic manipulation due to the more strict requirements of solutions. 
Due to these factors, the applications of symbolic and algebraic manipulation to 
engineering problems has became more popular recently. Some of the publications, such 
as those written by Madson, Smith and Hoff [16], Levi [17], Wilkins [18], Noor and 
Anderson [19], Korncoff and Fenven [20], Steinberg and Roache [21], will be 
discussed in more detail in the next paragraph. 

4. Application in the other fields — In addition to science and engineering, Symbolic and 
algebraic manipulation has been applicable in other fields, such as information 
management [22], education [23] [24] and business [25]. 

As time goes on, more and more applications will be reported in various fields. This is 
due to the fact that : 

1. The ongoing improvement in the memory space of hardware systems, especially 

personal computers. 

2. The availability of variously sound SAM software systems. 

However, in order to see that scratch paper is replaced by the computer screen in all 
areas, the people in the educational field should assume the responsibility of utilizing this new 
tool. As the discussion in the paper, written by Richard Pavelle in 1985, points out [25], only 
about 20 percent of people in the related field are aware of the existence of SAM system and 
less than a quarter of these actually use them. 

2.2 Reviews of SAM applications in engineering 

( 1) Computer algorithms for solving non-linear problems 

A paper [16] published in 1965 is believed to be the earliest document which employed 
computerized symbolic and algebraic manipulation to solve an engineering problem. The 
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authors, W. A. Madson, L. B. Smith and N. J. Hoff, developed their own software to find the 
solution for the post-buckling behavior of thin-walled circular cylindrical shells under axial 
compression. The major commands developed by them were : 

SERIESMULT . To expand the expressions, e.g. (a*sin(x)+b*cos(y)) n . 

*TRIGSPAND : To treat non-double trigonometric terms, such as sin 2 (x ), 
sin(x)*sirt(2x)*cos(y) etc., into double trigonometric terms like cos(y )*cos(x ). 

*SEARCHSTORE : To search and collect the coefficients of like trigonometric function, 
then store them. 


*NEWTNRAPH : To solve the nonlinear system equations obtained from the calling of the 
last three commands by using the Newton-Raphson iteration method. 

The application of the above commands to the shell post-buckling problem started at the 
assumption of radial displacement, 

w=t IA-jCos( inx/kx)cos(jjvy/Xy ) 

Then by the strain-displacement relationship and Hook's law, the stresses could be obtained. 
As the stresses (therefore the Airy stress function) were known, the membrane energy, 
bending strain energy and the potential of axial load could be derived. The resultant total 
potential energy was then minimized with respect to the coefficients of radial displacement w. 
The system equations obtained after the minimization then could be solved by calling the 
NEWTNRAPH command. 

(2) Symbolic algebra by computer-applications to structural mechanics 

One of the earliest publications of symbolic manipulation application in the engineering 
field was in 1971, when only a few SAM systems existed. Only REDUCE and FORMAC 
were mentioned in this paper. At the beginning of the paper [17] by 1. M. Levi, he described 
the story of SAM application in seeking the minimum theoretical post-buckling load for a thin 
circular cylindrical shell under axial compression. Starting from 1941, Von Karman and Tsien 
indicated that a low post-buckling load could be found with only two terms included in the 
series expression of normal displacement. This inconsistency in the Von Karman-Tsien's 
solution was not found until J. Kempner increased the series into three terms and found a 
further lowering of the post-buckling load in 1954. Since then, additional investigations were 
conducted as new terms were added into the series solution. But finally everybody was limited 
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by their inability to solve the complex algebraic equations without error. The drive to seek the 
minimum load did not end until 1965 when Madson, Smith and Hoff of Stanford University 
wrote the special program in ALGOL to increase the series into 14 terms and found the 
minimum load approached zero, which revealed the basic fallacies in the application of the 
Karman-Tsien procedure. 

The second topic in the paper talked briefly about the derivation of a stiffness matrix for 
a compatible triangular plate bending element by symbolic and algebraic manipulation. Then an 
example of the calculation of creep strain rate in plate and shell problems was demonstrated 
using SAM. This computation started from the x, y components of stress which were 
expressed as a double trigonometric series, followed by the calculation of equivalent stress, 
and ended in the substitution of the above quantities into strain rate equations. The resultant 
fortran codes were printed out by REDUCE. The paper ended with a brief discussion on the 
REDUCE capacities. 

(3) Applications of symbolic algebra manipulation language for composite structures analysis 

This paper [18] was published in 1973 by Dick J. Wilkins, Jr. of General 
Dynamics/Convair Aerospace Division, Fort Worth, Texas. The author used PL/I FORMAC to 
calculate the strain energy for an anisotropic shell. The strain energy can be written as 

", 

M x ’ 

My 

Where 

• N : stress resultants. 

• M : moment resultants. 

• e : mid-plane strain. 

• k : curvature. 

• £2: shell surface area. 
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• V : strain energy. 


and the constitutive equations are expressed as follows : 
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Where : Ay, By, Dy are the pertinent constitutive components of Hooke's Law. The 
equation (2. 1) and (2.2) can be combined into the form 

V ~ t^H £ } + 2{e} [£]{*•} + {*-} r [D]{K- (2.3) 


Then the displacements are approximated as 
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Where the Cy are constants to be determined by Rayleigh-Ritz method. By Vlasov shell 
theory, the strain-displacement relations are 
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Where R here is the radius of shell curvature. 


With the assumption of a symmetric constitutive matrix, the calculation of the integrand 
of equation (2.3) was done by FORMAC by the substitution of equations (2.4) to (2.12) into 
(2.3). The Rayieigh-Ritz method was then applied to take the partial derivative with respect to 
all undetermined constants in the displacement series. This was also manipulated by 
FORMAC. The resultant expressions were the energy variation which is in terms of ~ , — 
and their derivatives. Since FORMAC was incapable of performing the symbolic integration at 
that time, the informal "symbolic integration" was done by examining each term in the energy 
variation for a specific combinations of derivatives. Each time a certain type was found, it was 
replaced by a symbol and an appropriate constant to allow for the non-dimensionalization of the 
integrals. There were a total of twelve different integrations in the energy variation equation. 
The final expressions were then slightly modified into FORTRAN code by adding DO loops 
and suitably changing the indices by hand. 

(4) Computerized Symbolic Manipulation in Structural Mechanics — Progress and Potential 

In the beginning of the paper [19], A. K. Noor and C. M. Anderson introduced the 
symbolic and algebraic manipulator MACSYMA. These included the brief history, basic 
capacities and special commands, as well as associated packages. 

The second part of the paper gives three applications in the structural mechanics field by 
using MACSYMA. They are : 

1. Generation of characteristic arrays of finite elements for a shear flexible shallow shell 
element — There were three types of basic integrals for linear problems and three types 
of basic integrals for geometrically nonlinear problems. These were 

(a) Linear problems 

A'-f^N'N’dQ (2.13) 

Q 

B J ■ / (2.14) 

a 

c * -S a< .,d.N'd f N‘dQ (2.15) 

(b) Geometrically nonlinear problems 
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The evaluation of integrals in equations (2.13) and (2.14) can be performed 
analytically by MACSYMA. However in genetal the integrand of equattons (2. 15) to 
(2.18) cannot be integrated exactly due to the existence of a Jacobian determinant in 
the denominator of the integrand* Therefore the hybrid approach (numerical 
quadrature plus symbolic manipulation) was proposed. The number of integrations to 

be performed can be substantially reduced by the help of permutative and Dihedral 
symmetries. 

2. Evaluation of effective stiffness and mass coefficients of continuum models for repetitive 
lattice structures — The symbolic manipulations by MACSYMA included the evaluation 
of strain components, calculation of strain energy (with the thermoelastic strain energy) 
and kinetic energy, computation of stiffness and thermal coefficients as well as effective 
mass coefficients, forming the Lagrangian of the system and finally obtaining the 
governing differential equations. The numerical analysis started as soon as the governing 
equations were obtained. This numerical analysis was also done in MACSYMA. The 
results of mode shapes were then plotted out by MACSYMA 's graphic facility. 

3. Application of the Rayleigh-Ritz technique to the free vibration analysis of laminated 
composite elliptic plates — The tasks done by MACSYMA in this application were 

(a) Selecting approximation functions for each of the fundamental unknowns 
displacement amplitude with undetermined coefficients and developing analytic 
expressions for the specific strain and kinetic energies as quadratic functions of the 
undetermined coefficients. 

(b) Differentiating specific strain and kinetic energies with respect to the undetermined 
coefficients symbolically. 

(c) Evaluating stiffness and mass coefficients by performing integrations over volume. 


"This has been done successfully, see the details in chapter three. 



(d) Simplifying the expressions for the nonzero stiffness and mass coefficients and 
developing FORTRAN code. 

After stiffness and mass coefficients had been evaluated, the vibration frequencies and 
mode shapes could be obtained numerically by using any scheme for generalized eigenvalue 
problems. 

The last part of the paper discussed the problems which limited the applicability of 
computerized symbolic manipulation. The major problems mentioned in the paper were 
summarized as follows. 

1. Production of large expressions during the computation (intermediate expressions 
swell). 

2. Slow speed of symbolic computation. 

3. Low portability of large symbolic manipulation systems. 

4. Need for analyst interaction during the symbolic computation. 

5. Inability to estimate the storage requirements and CPU time for symbolic computations. 

6. Problems associated with interface between algebraic and numerical calculations. 

In addition, the authors suggested the directions of future research in this field. They 

were 

1. Reduction of a general (tensor) formulation of structural mechanics problem to its 
computational level. 

2. Hybrid computations. 

3. Approximate symbolic integration of rational functions. 

(5) Symbolic generation of finite element stiffness matrices 

As the title implies [20], the authors A. R. Komcoff (Boeing computer service, Seattle 
WA) and S. J. Fenves (Camegie-Mellon University) used the symbolic processor MACSYMA 
to assist in the development of a software to generate the stiffness matrices for finite element 
analysis. These included the construction of the strain-displacement matrix, calculation of the 
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determinant of the Jacobian, and multiplication of relevant matrices. The integrand was then 
integrated symbolically if it was integrate (e.g. the constant strain triangle element). Otherwise 
it would be output as the function of the problem parameters for further numerical evaluation 
(e. g. four-node quadrilateral element^ . In addition, the software gave two options in the 
material property matrix. One was "user-supplied" material property. The other was "library 
supplied" which provides one-dimensional elasticity, plain stress, plain strain, axisymmetric, 
and 3-D linear isotropic elasticity The example of constructing the isoparametric formulation 
for constant strain triangle was also shown in the appendix of the paper. 

(6) Symbolic manipulation and computational fluid dynamics 

The authors S. Steinberg and P. j. Roache [21] employed the symbolic and algebraic 
manipulator VAXIMA, a VAX version of MACSYMA, to transform the physical differential 
equation and the boundary conditions into the rectangular region and then constructed the so 
called stencil coefficient matrix for a finite difference scheme. The major ideas came from the 

general elliptic problems. In physical coordinates, the linear elliptic equation could be 
expressed as 
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Where a-, b- v c, d were given and were the function of coordinates in general. The problem 
was to find a numerical approximation solution which satisfies equation (2.20) and the given 
boundary conditions. 


Since the physical domain is not regular in general, it is necessary to transform the 
physical coordinates into the rectangular, computational coordinates in which the finite 
difference scheme could be constructed easily. This coordinate transformation involved the 
calculations of the Jacobian matrix, its determinant, and cofactors. The equation in new 
coordinates would become 


Cf- 



( 2 . 21 ) 


The integration (2.15) fora four-node isoparametrical quadrilateral element has been 
obtained exactly and will be discussed in detail in chapter three of this report. 
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Where the tilde denoted that the quantities were functions of computational coordinates. 


After the transformation of equation and boundary conditions had been done, the 
centered difference method was employed to construct the finite difference scheme as follows : 
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Where the coefficient c-j ^ were called the stencil and was constructed by symbolic 
manipulation. Taking advantage of the symmetric property, the number to be computed for 
c i,j,k cou ld dropped to 10 from 27. The resultant expressions of c— ^ then could be coded 
in the FORTRAN language for the next numerical scheme. 


2.3 Conclusion 


After making a survey of the publications on symbolic and algebraic manipulation, the 
following conclusions are drawn : 

1. None of the papers applying SAM to engineering problems tried to get closed-form 
solutions. They kept traditional methodology by increasing the terms of the 
approximation function to get more accurate solutions. This is due to the difficulty in 
solving generally partial differential equations or integral equations analytically. 

2. The papers discussing the application of symbolic and algebraic manipulation on the 
finite element analysis stop at the step of making a local stiffness matrix, local mass 
matrix, etc. The same situation also occurred in finite difference analysis. This was 
because 

(a) The finite element and finite difference methods are themselves approximation 
methods. The accuracy of results depends on many factors, not just on round-off 
error or integration error which can be cured by symbolic and algebraic 
manipulation. Although it was also one of the purposes to improve the accuracy of 
the solution, the major consideration in applying symbolic and algebraic 
manipulation was to help in the formulation of the tedious mathematical equations. 

(b) In general engineering problems, the stiffness matrix in FEA and stencil coefficient 
in FDA are huge in dimension. The limitation of memory space makes the execution 
of FEA's (or FDA) job impossible by symbolic and algebraic manipulation. 
Therefore it is necessary to be finished by numerical analysis. 
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(c) Although most of the SAM systems also possess the capacity of numerical analysis, 
the execution speed of numerical analysis in symbolic and algebraic manipulator is 
slower in comparison to that in pure numerical analysis. The difference of 
efficiencies between them is remarkable when the job is big. Therefore it is best not 
to have it done completely in symbolic and algebraic manipulation. As the documents 
showed, nobody did the whole FEA or FDA job in symbolic mode alone. 
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CHAPTER III 


CAPABILITIES OF THE SYMBOLIC AND ALGEBRAIC MANIPULATORS 


3.1 Introduction 

Rou^r T bi "“ eS ° f ' he Symb °" C ^ a ' 8ebraiC ma " ipUla ' 0r are 1 ulte s ^"> ^"dent. 

Roughly speaking, a system which is destgned for general purpose usage usually possesses 
e fune, tons of different, atton, integration, mainx operation, polynomial manipulation, pattern 

MA™ “' ,Ut, 7 an f d h eq,,ati0n SO ' Ver ' S ° me «■> “ MACSYMA and 

EMATICA, have a lot of built-in mathematical functions which allow the users to get 

answers by just calling the appropriate command once. Others, like REDUCE may need 
users to wn te a short program to get the same answers. 

This chapter will demonstrate the fundamental capabilities of the symbolic and algebraic 

^pulators Which are available a, hand by solving examples of applied mechanics. Since 

REDUCE .S the oldest system available at The University of Michigan, most examples will be 

demonstrated by using REDUCE. Of course, MACSYMA w,ll be employed to help the 
demonstration if it is necessary. 

MACSY U MA r r a,ely ' r UCE d0eS "'‘ POSSeSS the graphiC «» version of 

ACSYMA being used ,n The University of Michigan also doesn't include the graphics 

package, although ,t ,s available on the market. Therefore, the postprocessing of the results 

from symbolic and algebraic manipulators w,ll be done by other graphics packages. 

3.2 What can the symbolic and algebraic manipulators do ? 

In this section, some of the most useful operations in symbolic and algebraic 
manipulation are demonstrated in detail by examples. They are differentiation, integration 
matrix operation, algebraic equation solving, treatment of trigonometric fund, on, different! 
equation solving, polynomial and rational operation, fortran code output, number system 

substitution and built-in functions. The strategies and parttcular techniques are also mentioned 
at the place where they are necessary. 
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3.2. 1 Differentiation 


All differentiations, without exception, can be done analytically by REDUCE If it is 
necessary, REDUCE knows how to apply the chain rule to solve problems. The powerful 
capabiht.es of this analytical differentiation will probably replace traditional numerical 
differentiation in many cases, such as the evaluations of the Jacobian and Hessian matrices 

However, the wrong results may be obtained by careless or naive users. For instance, ,n 
finding the first derivative of f wiIh respec , lo x Two dlffere „, ^ obtaj[]ed ^ 


1: on time; 

Time: 134 ms 

2: df(x**(x**x),x); 

X 

X 2 

X +X*(LOG(X) *X + LOG(X)*X + 1) 
X 

Time: 383 ms 

3: df(x**x**x,x); 

2 

X 

X *X*(2*LOG(X) + 1) 

Time: 233 ms 


Here the first solution is correct. How to judge the correctness of the results is one of 
the important tasks in symbolic manipulation. A sound background knowledge in the SAM and 
problem-related fields is very helpful in checking them. 

In some cases, the unevaluated differentiation form of functions, such as is desired 
to be retained throughout the computation. This also can be done as follows : 

4: depend x,t; 

Time: 84 ms 

5: depend y,t; 

Time: 83 ms 

6: p:=a*x*y; 

P := A*X*Y 
Time: 150 ms 

7: df(p,t); 

A*(DF(X,T)*Y + DF(Y,T)*X) 

Time: 133 ms 
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3.2.2 Integration 


In REDUCE, all the integrations performed by the command INT are indefinite 
integrations with the integration constants discarded. If the function is not integrable by 
REDUCE (the closed-form solution may exist theoretically), the original form will be displayed 
on the screen. The definite integration may be obtained by further substituting the upper and 
lower limits into the results after the indefinite integration. 

%A non-integrable case. 

8: int(sqrt(a A 2-x A 2),x); 

2 2 

INT(SQRT(A - X ),X) 

Time: 950 ms 

%An integrable case. 

9: int(l/(a A 2+x A 2),x); 

X 

ATAN( — ) 

A 


A 

Time: 466 ms 

In most cases, the non-integrable integrand will become integrable after appropriate 
manipulation. This pre-treatment involves the technique of changing the integrating variables in 
fundamental calculus. Sometimes the intelligent users can substantially extend the capabilities 
of the symbolic and algebraic manipulator by suitably combining human intelligence with the 
tireless and errorless advantages of computer. For example, if the x in command 8 is 
substituted by a*cos(t), then dx=-a*sin(t)*dl and the integration of V « 2 -* 2 with respect to x 
will become the integration of -a?*sir?(t)dtmih respect to/, which is integrable by REDUCE. 
After the integration is done, the original variable jc may be substituted back to get the desired 
expression in terms of x. The check may be done by skeptics by differentiating the resultant 
expression to get the original integrand. The following three commands demonstrate these 
procedures. 

10: int(-a A 2*(sin(t)) A 2,t); 

2 

A *(COS(T)*SIN(T) - T) 

2 

Time: 2134 ms 
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11 : sub(cos(t)=x/a,sin(t)=sqrt(l-x** 2 /a** 2 ),t=acos(x/a),ws); 

X 2 2 2 

ACOS(— )*A - SQRT(A - X )*X 
A 


2 

Time: 716 ms 

12: df(ws,x); 

2 2 
A -X 


2 2 
SQRT(A -X) 
Time: 450 ms 


3.2.3 Matrix operation 

Matrix operation is one of the powerful capabilities of symbolic and algebraic 
manipulators. These include the addition and multiplication of matrices, multiplication of 
matrices and scalars, inverting matrices, calculating the determinant of a square matrix, finding 

the trace, computing the eigenvalues and associated eigenvectors exactly if they are available 
and so forth. 


3.2.3. 1 Matrix multiplication 

The three body ngtd rotation 1-2-3 in dynamics is a good example of the utility of 
matnx multiplication. In robotics, it is necessary to find the analytical foim of final orientation 
from which the rotation angles of each arm can be computed. The final direction cosine matrix 
d is obtained from the product of three consecutive direction cosine matrices a, b, c. 

%Declaring four matrices. 

13: matrix a(3,3),b(3,3),c(3,3),d(3,3); 

%Inputting matrices. 

In a i : \ = ^ a i t((1,0,0) ^ 0,COS(ql), ' sin(ql)),(0,sin(< l 1 )* cos ^ 1 )))i 

A( 1,1) 1 

A(l,2) := 0 
A(l,3) := 0 
A(2,l) := 0 
A (2,2) := COS(Ql) 

A(2,3) :=-SIN(Ql) 

A(3,l) := 0 
A(3,2) := SIN(Ql) 

A(3,3) := COS(Ql) 
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Time: 700 ms 

15: b:=mat((cos(q2),0,sin(q2)),(0,l,0),(-sin(q2),0,cos(q2)))$ 

Time: 367 ms 

16: c:=mat((cos(q3),sin(q3),0),(-sin(q3),cos(q3),0) ,(0,0,1))$ 

Time: 383 ms 

%multiplication of three matrices. 

17: d:=a*b*c; 

D(l,l) :=COS(Q2)*COS(Q3) 

D( 1 ,2) :=COS(Q2)*SIN(Q3) 

D(l,3) :=SIN(Q2) 

D(2,l) :=-(COS(Ql)*SIN(Q3)-COS(Q3)*SIN(Ql)*SIN(Q2)) 
D(2,2):=COS(Ql)*COS(Q3)+SIN(Ql)*SIN(Q2)*SIN(Q3) 

D(2,3) :=-COS(Q2)*SIN(Ql) 

D(3,l) :=-(COS(Ql)*COS(Q3)*SIN(Q2)+SIN(Ql)*SIN(Q3)) 
D(3,2) :=-(COS(Ql)*SIN(Q2)*SIN(Q3)-COS(Q3)*SIN(Ql)) 
D(33) :=COS(Ql)*COS(Q2) 

Time: 617 ms 


The terminators $ in command lines 15 and 16 prohibit the printing of results and save 
almost half of the time compared to command line 14 which uses the other terminator. 

3. 2.3. 2 Matrix inversion 

Unlike numerical analysis in which the time-consuming operation of matrix inversion is 
to be avoided, to find the inverse of a matrix symbolically is one of the significant and simple 
tasks in symbolic and algebraic manipulation. One important application of it is in solving a 
system of linear equations. For example, the problem of finding a curve to fit the given set of 
data by the least square method results in solving a system of linear equations. The coefficient 
matrix here is the Hilbert matrix which is usually used to investigate the phenomenon of round- 
off error accumulation. REDUCE can solve this problem exactly. The numerical solution and 
symbolic solution are tabulated in Table 3.1, 3.2, 3.3 for three different Hilbert matrix sizes. 
As the tables show, the numerical solution is not capable of producing accurate results even for 
the case of the 7*7 Hilbert matrix. The deviation between both solutions is also plotted in 
Figure 3.1. The significance of symbolic and algebraic manipulation is evident. 

8: matrix h(40,40),x(40,l)$ 

9: for i:=l:40 do forj:=l:40 do h(i,i):=l/(i+j-l)$ 

Time: 30917 ms 

10: for i:=l:40 do x(i,l):=i$ 

Time: 583 ms 

11: h:=(l/h); 
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H(l,l) :=1600 
H( 1 ,2) :=- 1279200 
H(l,3) : =340267200 
H(l,4) : =-45 113759600 
H( 1 ,5) : =3573009760320 
H( 1 ,6) :=- 187583012416800 


H(40,39) :=-l 141 149866470104951399125616277120066810096080000 
H(40,40) := 58520505972825894943544903398826670092825440000 
Time: 1355067 ms 

12: x:=h*x; 

X(l,l) := -64000 

X(2,l) := 102272040 

X(3,l) := -40781023920 

X(4,l) := 7204667408120 

X(5,l) := -712815447183840 

X(6,l) := 44879235720719400 

X(7,l) := -1948531047013671120 

X(8,l) := 61638279321584754360 

X(9,l) := -1478390066741437724160 

X( 10, 1 ) : = 2770696 1 874232024704320 

X(1 1,1) := -415343205971360018352000 

X(12,l) := 5073605405648309638180800 

X(13,l) := -51267503668234803803526400 

X(14,l) := 433831946060133824442926400 

X(15,l) := -3105695134246771063279900800 


X(33,l) := -275148980879869194392659362117120 
X(34,l) := 129027970745724768036578024940720 
X(35,l) := -49525830202172298308155472545440 
X(36,l) := 15151287095628987186696245706000 
X(37,l) := -3551734684918636715496119886400 
X(38,l) := 598923394712817307441549441200 
X(39, 1) := -64662238360272798660726511200 
X(40,l) := 3356375056654755429131776400 
Time: 96000 ms 
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X 

Exact REDUCE solution 

Gauss elimination solution 

X(l) 

125 

124.7899166408321 

x(2) 

-2880 

-2875.997324887610 

x(3) 

14490 

14472.49323423709 

x(4) 

-24640 

-24613.29019994230 

x(5) 

13230 

13216.83350607823 

Tal 

3le 3. 1 : Comparison of solution for 5*5 Hilbert matrix 

X 

Exact REDUCE solution 

Gauss elimination solution 

x(l) 

-216 

-204.4675087167038 

x(2) 

7350 

7027.565254454758 

x(3) 

-57120 

-54968.63675793860 

x(4) 

166320 

160780.3224850843 

x(5) 

-201600 

-195532.2818417542 

x(6) 

85932 

83555.64156991533 

Tat 

>le 3.2 : Comparison of solution for 6*6 Hilbert matrix 

X 

Exact REDUCE solution 

Gauss elimination solution 

x(l) 

343 

131.3201346851223 

x(2) 

-16128 

-7941.234641235595 

x(3) 

177660 

100320.8278414759 

x(4) 

-772800 

-476154.4030660979 

x(5) 

1559250 

1020735.514691367 

x(6) 

-1463616 

-1001760.776649569 

x(7) 

516516 

365761.5297697352 


Table 3.3 : Comparison of solution for 7*7 Hilbert matrix 
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Figure 3.1 : Error between the solutions of REDUCE and Gauss elimination 


3.3 Eigenvalues and Eigenvectors 

To find the eigenvalues and eigenvectors for a matrix is another important task in the 
application of symbolic and algebraic manipulation. Since to solve the eigenvalue problems 
analytically for an arbitrary dimensional matrix is theoretically impossible, the following 
examples only show the exact solutions for a 3 by 3 matrix. 

21: matrix s(3,3)$ 

22: s:=mat((sxx,sxy,sxz),(sxy,syy,syz),(sxz,sxy,szz))$ 

23: mateigen(s,eta); 

3 2 2 2 ? 

{{ETA-ETA *SXX-ETA *SYY-ETA *SZZ+ETA*SXX*SYY+ETA*SXX*SZZ-ETA*SXY~ 

2 

-ETA*SXY*SYZ-ETA*SXZ +ETA*SYY*SZZ+SXX*SXY*SYZ-SXX*SYY*SZZ 
2 2 2 
-SXY *SXZ+SXY *SZZ-SXY*SXZ*SYZ+SXZ *SYY, 

1 , 
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ARBC0MPLEX(1)*(ETA*SXZ+SXY*SYZ-SXZ*SYY) 

MAT(1,1):= 

2 2 

ETA -ETA*SXX-ETA*SYY+SXX*SYY-SXY 

ARBC0MPLEX(1)*(ETA*SYZ-SXX*SYZ+SXY*SXZ) 

MAT(2,1):= ' 

2 2 

ETA -ETA*SXX-ETA*SYY+SXX*SYY-SXY 

MAT (3 , 1 ) : = A RBCOMPLEX( 1 ) } > 

Time: 1283 ms 

24: s:=mat((5,l,0),(l,2,4),(0,4,3))$ 

Time: 284 ms 

25: mateigen(s,eta); 

3 2 

{{ETA -10* ETA + 14*ETA+53, 

1, 

4*ARBCOMPLEX(2) 

MAT(1,1) := 

2 

ETA - 7*ETA + 9 

4*ARBCOMPLEX(2)*(ETA - 5) 

MAT(2,1) := — — — 

2 

ETA - 7*ETA + 9 

MAT(3,1) := ARBCOMPLEX(2)}} 

Time: 666 ms 

26: trace(s); 

SXX + SYY + SZZ 
Time: 100 ms 


As the command lines 23 and 25 show, the solutions from REDUCE by calling 
MATEIGEN contain three parts. They are 


(a) Characteristic equation. 


(b) The number of repetition roots. It's one in the above examples. 


(c) Eigenvectors. 
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If the eigenvalues are required, the characteristic equation needs to be solved in 
addition. The equation solver command SOLVE which will be demonstrated later can meet this 
requirement. The ARBCOMPLEX(l) which appeared in the eigenvectors is referred to as 
"arbitrary complex constant". Examples 24 and 25 also show results from REDUCE by 
substituting numbers into matrix S. The trace of the matrix is also evaluated by the command 
TRACE as shown in example 26. 

3.2.4 Equation solver 

REDUCE can solve the polynomial equation up to order three exactly. If the solution 
includes the imaginary part, the "I" will show up to represent the imaginary symbol. The 
following example solves the characteristic equation obtained above. 

27: solve(ETA**3-10*ETA**2+14*ETA+53=0,eta); 

2/3 

{ ETA =-( (63 * SQRT ( 229) * I -69 1 * SQRT (3 )) *SQRT(3)*I+(63*SQRT(229)*I-691* 

„ 2/3 1/3 1/3 1/6 2/3 

SQRT(3)) -20*(63*SQRT(229)*I-691*SQRT(3)) *2 *3 -58*2 *SQRT( 3 )* 

1/3 2/3 1/3 1/3 1/3 i /g 

*3 *1+58*2 *3 )/(6*(63*SQRT(229)*I - 691*SQRT(3)) *2 *3 ), 

2/3 

ETA=((63*SQRT(229)*I - 691 *SQRT(3)) *SQRT(3)*I - (63*SQRT(229)*I - 691* 

2/3 1/3 1/3 1/6 2/3 

SQRT(3)) +20*(63*SQRT(229)*I-691*SQRT(3)) *2 *3 -58*2 *SQRT( 3 )* 

1/3 2/3 1/3 1/3 1/3 1/6 

3 *1 - 58*2 *3 )/(6*(63*SQRT(229)*I - 691*SQRT(3)) *2 *3 ), 

2/3 1/3 

ETA=((63*SQRT(229)*I - 691*SQRT(3)) + 10*(63*SQRT(229)*I - 691*SQRT(3)) 

1/3 1/6 2/3 1/3 1/3 1/3 1/5 

*2 *3 +58*2 *3 )/(3*(63*SQRT(229)*I - 691*SQRT(3)) *2 *3 )\ 

Time: 5717 ms 

If the numerical mode NUMVAL, complex switch COMPLEX and FLOAT mode are 
turned on, the numerical solution of three eigenvalues can be obtained in sixteen digits 
precision by default. The imaginary parts in the following example are very small and are due 
to the round-off errors. 

28: solve(ETA**3 - 10*ETA**2 + 14*ETA + 53=0,eta); 

{ETA=4.83075995061 1553d0 + 2.991 12422333 1416d-7*I, 
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ETA=-( 1.6167630306368408d0 + 6.960060167347475d-8*I), 
ETA=6.786003556862447d0 + (-2.2951 18206596669d-7)*I> 
Time: 1567 ms 


REDUCE can solve systems of linear algebraic equations exactly. The limitation is 
determined only by the memory capacity of the hardware system. The following example finds 
the minimization of a quadratic function 


Q=kl*yl*{2}+k2*yl*y2+k3*y2*{2}+k4*y2*y3+k5*y3*{2}-k6*y3 subject to yl+y2=2 

29: Q:=kl*yl**2+k2*yl*y2+k3*y2**2+k4*y2*y3+k5*y3**2-k6*y3+y4*(yl+y2-2)$ 
Time: 550 ms 

30: a:=df(q,yl); 

A := 2*K1*Y1 + K2*Y2 + Y4 
Time: 167 ms 

31: b:=df(q,y2); 

B := K2*Y 1 + 2*K3*Y2 + K4*Y3 + Y4 
Time: 167 ms 

32: c:=df(q,y3); 

C := K4*Y2 + 2*K5*Y3 - K6 
Time: 166 ms 

33: d:=df(q,y4); 

D := Y1 + Y2 - 2 
Time: 150 ms 


34: solve({a=0,b=0,c=0,d=0},{yl,y2,y3,y4}); 


2 

4*K2*K5 - 8*K3*K5 + 2*K4 - K4*K6 

{{ Y1= 

2 

4*K1*K5 - 4*K2*K5 + 4*K3*K5 - K4 
8*K1*K5 - 4*K2*K5 - K4*K6 


> 

2 

4*K1*K5 - 4*K2*K5 + 4*K3*K5 - K4 


Y3= - 


2*(2*K1*K4 - K1*K6 - K2*K4 + K2*K6 - K3*K6) 


2 

4*K1*K5 - 4*K2*K5 + 4*K3*K5 - K4 


2 2 

16*K1*K3*K5-4*K1*K4 +2*K1*K4*K6-4*K2 *K5-K2*K4*K6 


2 

4*K1*K5 - 4*K2*K5 + 4*K3*K5 - K4 


Y4= - 


}} 



Time: 1400 ms 


After the substitution of yl to y4 into the quadratic function, the minimum is found as 
follows : 


35: q:-q; 

2 2 2 

Q:=(16*K1*K3*K5-4*K1*K4 +4*K1*K4*K6 -K1*K6 -4*K2 *K5-2*K2*K4*K6+ 

2 2 2 
K2*K6 -K3*K6 )/(4*Kl*K5-4*K2*K5+4*K3*K5-K4 ) 

Time: 267 ms ’ 

3.2.5 Treatment of the trigonometric function 

REDUCE doesn't even know an equation as simple as sin 2 (q)+cos 2 (q)=l. However 
REDUCE does possess the potential to learn it. Due to this powerful capability, the 
trigonometric functions can be handled easily by just teaching REDUCE the operation rules. 
For example, without teaching the operation rule of trigonometry, the determinant of direction 
cosine matrix d in command 17 is as follows : 

36: det(d); 

2 2 2 2 2 2 2 
COS(QU *COS(Q2j *COS(Q3) +COS(Ql) *COS(Q2) *SIN(Q3) +COS(Ql) * 

COS(Q3) *SIN(Q2)^+COS(Ql) *SIN(Q2) *SIN(Q3) _ +COS(Q2) *COS(Q3)i * 

SIN(Q1HC0S(Q2)^*SIN(Q1) *SIN(Q3) +COS(Q3) *SIN(Q1) 2 *SIN(Q2) 2 + 

SIN(Q1)“*SIN(Q2) *SIN(Q3) 

Time: 483 ms 

After teaching REDUCE the appropriate operation rules, the solution becomes quite 
simple. Note that the time consumption^ command 38 is longer than that in command 36. 

This is due to the extra work needed for simplification. It also reveals the phenomenon of 
internal swells. 

Tim^^W ms^** 2+S * n(ql) ** 2 ~ 1,COS(q2) * ll!2+Sin(q ^ ) ** 2=1,CC>S(q3) ** 2+S * n ^^**^ = ^’ 

38: det(d); 

1 

Time: 584 ms 
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3 . 2.6 Solving differential equation 


REDUCE is unable to solve the differential equation directly, while MACSYMA does 
possess this capability. The following example is a problem of beam deflection w(x) under 
uniform load q [Figure 3.2], The differential equation is of the form 


d 2 vv 
d jc 2 


- 5 * w + r*x (x - L) 


(3. 1) 


Where 5 and r ,in general, are the function of Young's modules, moment of inertial as 
well as boundary conditions. In the case of small deflection with simple supported on both 
end, the s becomes zero, and r=q/2EI. 

w . 

q 


rrrrnrrrrrrrn 



^ l H 


Figure 3.2 : Beam under uniform load. The boundaries are not 
specified to find general solution. 

The (Cn) in the following examples is the MACSYMA prompt for inputting the 
command and (Dn) is the solution given by MACSYMA. 

(Cl) depends(w.x); 

(Dl) [W(X)] 

(C2) diff(w,x,2)-s*w-r*x*(x-l)=0; 

2 

d W 

(D2) - R X (X - L) + S W = 0 

2 

dX 

(C3) ode2(d2,w,x); 

Is S positive, negative, or zero? 

p; 
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2 

SQRT(S) X -SQRT(S) X RSX-LRSX+2R 

(D3)W = %K1 %E + %K2 %E 

2 

S 


Where the %K1 , %K2 are constants to be determined by boundary conditions. The 
%E here is the symbol of exponential function. 


Although the next example (no specific physical problem associated with it) is more 
complex, it only takes 3.6 milliseconds in MACSYMA. 

(C4) depends(y,x); 

(D4) [Y(X)] 

(C5) diff(y,x,2)+(2/x)*diff(y,x)-(2/x**2)*y-(l/x**2)*sin(log(x))=0; 

dY 

2 2 — 

d Y dX 2Y SIN(LOG(X)) 

(D5) + = 0 

2 X 2 2 

dX XX 

(C6) ode2(d2,y,x); 

3 SIN(LOG(X)) + COS(LOG(X)) %K2 

(D6) Y = + %K1 X + 

10 2 

X 

(C7) time(d3); 

Time: 

(D7) [3.6d0] 


3.2.7 Polynomial and rational operations 


Polynomial and rational operations is one of the most important and useful functions in 
symbolic and algebraic manipulation. There are two occasions in employing these functions. 
First, in most cases the problems to be solved are not as simple as the above demonstrations. 
Therefore it is necessary to manipulate the formulae into a machine manageable forms before 
calling the appropriate REDUCE commands to solve them. Second, sometimes the solutions 
are restricted to specific forms for particular usage. In order to get the appropriate forms, there 
is no way to avoid employing this package. The capabilities of this package include the controls 
of the expansion, factorization, and cancellation of common factors, determining the GCD of 
two polynomials, obtaining the part of polynomial and rational functions, and so forth. 
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%Tuming off expansion switch. 
39: off exp; 


%Inputting p and q polynomials. The numerical common factors 
%will be automatically factored out. 

40: p:=(x-l)*(5*x-3)**2*(4*x+8) A 3*(9*x-6); 

2 3 

P := 192*(5*X - 3) *(3*X - 2)*(X + 2) *(X - 1) 

41: q:=(x-l)*(5*x-3)*(4*x+8)*(x+4); 

Q := 4*(5*X - 3)*(X + 4)*(X + 2)*(X - 1) 

%Getting greatest common divider of p and q. 

42: gcd(p,q); 

4*(5*X - 3)*(X + 2)*(X - 1) 

%Tuming on expansion switch and checking p, q. 

43: on exp; 

44: p; 

7 6 5 4 3 2 

192*(75*X +235*X -163*X -723*X +392*X +664*X -624*X+144) 

45: q; 

4 3 2 

4*(5*X + 22* X - 5*X - 46* X + 24) 

%Defining a fraction r. 

46: n=p/q; 

7 6 5 4 3 2 

R :=(48*(75*X +235*X -163*X -723*X +392*X +664*X -624*X+144))/ 

4 3 2 

(5*X +22*X -5*X -46*X+24) 

%Tuming on the greatest common divider switch and rechecking r. 

%The common factors have been cancelled as shown in command 47. 

47: on gcd; 

48: r; 

4 3 2 

48*(15*X + 41*X - 10*X - 52*X + 24) 


X + 4 

^Getting the denominator and numerator of fraction r. 
49: den(r); 

X + 4 
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50: num(r); 


4 3 2 

48*(15*X + 41*X - 10*X - 52*X + 24) 

%Getting the leading degree of numerator of fraction r 
51: deg(num(r),\); 

4 

S T on Sn‘° r ° f r in COmm “ d 49 a,SO ' “ factt)red “ «* multiplicaUons of each facto 

53: factorize(num(r)); 

{2,2,2,2,3,X + 2,X + 2,3*X - 2,5*X - 3} 

3.2.8 Fortran code output 


REDUCE can automatically produce the fortran expression, natural style expression 
(default), and REDUCE code. The fortran code can be made as a subroutine and be directly 
input to the fortran main program. The natural style expression allows it to be looked as hand- 
written form, while REDUCE code is useful in making a REDUCE subroutine for input into 
the REDUCE main program. Since the results from REDUCE are generally very lengthy, the 
functions of code-conversion make the switch from symbolic and algebraic manipulation to 
numerical analysis smoother. This not only saves effort in symbolic and algebraic 
manipulation, but also rules out all the possibilities of error introduced by hand typing. The 
following examples will give a clearer understanding about these functions. 

= t ( t a+b-c) e *^ lyn0mial ’ ThC ° UtpUt f0miS 316 * n natUral Sty,e of human writing. 


P:- A +7* A *B-7*A *C+21*A *B -42* A *B*C+21*A 5 *C +35*A *B 3 -105*A 4 *B 2 *C+105* 

4 243 34 33 322 3 T * a 

A *B*C -35* A *C +35*A *B -140*A *B *C +210*A *B *C -140*A *B*C +35*A *C 


2 5 2 4 2 3 

+21*A *B -105*A *B *C +210*A *B 


2 223 24 25 

*C -2 10* A *B *C +105* A *B*C -21* A *C +7* 


6 5 4 2 33 24 S * 

A*B -42*A*B *C+105*A*B *C -140*A*B *C +105*A*B *C -42*A*B*C +7*A*C 

7 6 52 43 34 25 67 

+B -7*B *C+21*B *C -35*B *C +35*B *C -21*B *C +7*B*C -C 
Time: 1434 ms 


%Tuming on the fortran-code conversion switch 
55: on fort; 
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%Checking the output of fortran code. 

56: p; 

ANS=A**7+7.*A**6*B-7.*A**6*C+21.*A**5*B**2-42 *A**5*B 
• *C+ 21.*A**5*C**2+35.*A**4*B**3-105.*A**4*B**2*C+105 
-*A**4*B*C**2-35.*A*M*C**3+35.*A**3*B**4-140 *A**3*B 
.**3*C+210.*A**3*B**2*C**2-140.*A**3*B*C**3+35.*A**3 
,*C**4+21.*A**2*B**5-105.*A**2*B*M*C+210.*A**2*B**3*C 
.**2-210.*A**2*B**2*C**3+105.*A**2*B*C**4-21 *A**2*C 
. **5+7.*A*B**6-42.*A*B**5*C+105.*A*B**4*C**2-140 *A*B 
**3*C**3+105.*A*B**2*C**4-42.*A*B*C**5+7 *A*C**6+B** 

. 7-7.*B**6*C+21.*B**5*C**2-35.*B*M*C**3+35 *B**3*C** 

. 4-21.*B**2*C**5+7.*B*C**6-C**7 
Time: 850 ms 

%Changing the number of continuation line in fortran code. 

57: cardno!*:=10$ 

58: p; 

ANS1=7.*A*B**6-42.*A*B**5*C+105.*A*B**4*C**2-140 *A*B 
. **3*C**3+105.*A*B**2*C**4-42.*A*B*C**5+7 *A*C**6+B** 
7-7.*B**6*C+21.*B**5*C**2-35.*B*M*C**3+35 *B**3*C** 

. 4-21.*B**2*C**5+7.*B*C**6-C**7 
ANS=A**7+7.*A**6*B-7 *A**6*C+21.*A**5*B**2-42 *A**5*B 
. *C+21.*A**5 !k C**2+35.*A*M*B**3-105.*A**4*B**2*C+105 
* A *M*B*C**2-35.*A*M*C**3+35.*A**3*B*M-140.^**3*B 
. **3*C+210.*A**3*B**2*C**2-140.*A**3*B*C**3+35 * A **3* 

. C**4+21.*A**2*B**5-105.*A**2*B*M*C+210.*A**2*B**3*C 
. **2-210.*A**2*B**2*C**3+105.*A**2*B*C*M-21 *A**2*C 
. **5+ANSl 
Time: 1034 ms 

%Tuming off fortran-code conversion switch. 

59: off fort; 

%Tuming off the 'natural style' function switch. 

60: off nat; 

%Checking the output of REDUCE code. 

61: p:=p; 

P := A**7+7*A**6*B-7*A** 

6*C+21*A**5*B**2-42*A**5*B*C+21 :,! A**5*C**2+35*A*M*B**3-105* 

A **4* B **2*C+105*A**4*B*C**2-35*A**4*C**3+35*A**3*B**4-140*A 

** 3 * B ** 3 * c+ 2 io*A** 3 *B** 2 *C** 2 - 140 *A** 3 *B*C** 3 + 35 *A** 3 *C** 4 + 

2i* a **2*b**5-105*A**2*B**4*C+210*A**2*B**3*C**2-210*A**2*B 

**2*C**3+105*A**2*B*C**4-21*A**2*C**5+7*A*B**6-42*A*B**5*C 

+105*A*B**4*C**2-140*A*B**3*C**3+105*A*B**2*C**4-42*A*B*C**5 

+7*A*C**6+B**7-7*B**6*C+21*B**5*C**2-35*B*M*C**3+35*B**3*C 

*M-21*B**2*C**5+7*B*C**6-C**7$ 

Time: 816 ms 
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3.2.9 Number system 


In addition to symbolic manipulation, the symbolic and algebraic manipulators can also 
do numerical analysis. There are three ways to treat numbers in REDUCE. They are 

(a) Integer — In general, there is no practical limit on the number of digits. For example, 
the value of 2 1000 gives 255 digits. It only takes 383 milliseconds. 

62: A:=2**1000; 

A : = 107 1 508607 18626732094842504906000181056140481 170553360744375038837 
0351051 124936122493 198378815695127594672917553 1468251871452856923 1 
4043 598457759857480393456777482423098542 1 07460506237 11418771821530 
464749835819412673987675591655439460770629145711 
Time: 383 ms 

(b) Fraction number — Numbers that aren't integers and operated with symbols (or 
numbers) are represented by default by the quotient of two integers with message(s) 
telling users the conversion. 

63: a:=0.999*b*c; 

*** 0.999 represented by 999/1000 
999*B*C 

A := 

1000 
Time: 200 ms 

64: a: =0.999*0.5; 

*** 0.999 represented by 999/1000 
*** 0.5 represented by 1/2 
999 

A := 

2000 

Time: 217 ms 


(c) Real number — It is also possible to ask REDUCE to work the floating point 
approximations to numbers with arbitrary precision with specified numbers of digit. 

%T uming on the numerical mode. 

65: on numval; 

%Tuming on the floating system switch. 

66: on float; 

67: pi; 

3. 141592653589793d0 
Time: 150 ms 

68: on bigfloat; 

*** Domain mode FLOAT changed to BIGFLOAT 
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% Specified 50 digits. 

69: precision 50$ 

70: pi; 

3. 141 59265 35897 93238 46264 33832 79502 88419 71693 993751 
Time: 217 ms 


3.2.10 Substitution 


There are two substitution functions in REDUCE. One is for local substitution, and the 
other is for global substitution. The difference can be revealed in the following examples. 

%Defining a function f. 

71: f:=6.*A**2*Rl**2*U2-3.*ATAN(Ul/A)*A**2* 

R1*R2*U1-9.*ATAN(U1/A)*A**2*R1*R2*U2+3.*ATAN(U1/A)*A 

**2*R2**2*U1+3.*ATAN(U1/A)*A**2*R2**2*U2-2.*ATAN(U1/ 

A)*R1**2*U1**3+6.*ATAN(U1/A)*R1**2*U1**2*U2$ 

Time: 1234 ms 

%making a local substitution and calling it as B. 

72: b:=sub(atan(ul/a)=k,0; 

2 2 2 2 2 
B :=-(3*A *K*R1*R2*U1+9*A *K*R1*R2*U2-3*A *K*R2 *U1-3*A *K* 

2 2 2 2 3 2 2 

R2 *U2 - 6* A *R1 *U2 + 2*K*R1 *U1 - 6*K*R1 *U1 *U2) 

Time: 433 ms 

%Checking the original function f after the local substitution. It's unchanged. 

73: f; 


U1 2 U1 2 U1 2 2 

-(3*ATAN( — )*A *R1*R2*U1+9*ATAN( — )*A *R1*R2*U2-3*ATAN(-— )*A *R2 *U1 
AAA 

U1 2 2 U1 2 3 U1 2 2 

-3*ATAN( — )*A *R2 *U2+2*ATAN( — )*R1 *U1 - 6*ATAN( — )*R1 *U1 *U2 
A A A 

2 2 

- 6* A *R1 *U2) 

Time: 367 ms 

%Making global substitution. 

74: letatan(ul/a)=g; 

Time: 133 ms 

%The original function f has been changed. 

75: f; 
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2 2 2 2 2 2 
'(3* A *G*R1*R2*U1+9*A *G*R1*R2*U2-3*A *G*R2 *U1-3*A *G*R2 *U2 

2 2 2 3 2 2 

- 6* A *R1 *U2 + 2*G*R1 *U1 - 6*G*R1 *U1 *U2) 

Time: 283 ms ; 


3.2.11 Built-in functions 


The built-in functions are quite system dependent. Since REDUCE is designed for 
general purpose usage, there are not many built-in functions. However they can be obtained by 
suitably combined commands of REDUCE. On the contrary, MACSYMA has many built-in 
functions which allow users to simply call commands once to get the solution. Some of these 
MACSYMA functions are shown in the following paragraphs. 

(a) Limit evaluation — If it is necessary, the function LIMIT in MACSYMA will 

automatically apply L'Hospital's rule to evaluate the formulae. 

(C8)limit(sin(x)/x,x,0,plus); 

(D8) K 1 

(C9) time(d8); 

Time: 

(D9) [2. 1 16d0] 

(CIO) limit((6*x*2+3*x-4)/(x-l),x,l,plus); 

(DIO) INF 

(Cl 1) limit((l-x)**(l/x),x,0); 

(DU) %E 1 

(C12) time(dll); 

Time: 

(D12) [2.75d0] 


(b) Laplace transformation - The LAPLACE command in MACSYMA can transform the 
functions in physical domain, such as EXP, LOG, SIN, COS, SINH, COSH, DELTA 
and ERF, into the j domain. In addition, it also can transform a differential equation 

into algebraic equation. The command for inverse of Laplace transform is also 
available. 


(C13) laplace( l/sqrt(t),t,s); 

SQRT(%PI) 

(Di3) 


SQRT(S) 
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(C 14) time(dl3); 

Time: 

(D14) [0.05d0] 

(C15)laplace(((c-b)*exp{a*t)+(a-c)*exp(b*t)+(b-c)*exp(c*t))/((a-b)*(b-c)*(c-a)),t,s); 

B-C A-C C-B 

+ + 

S-C S-B S-A 

(D15) 

(A - B) (B - C) (C - A) 


(C16) time(d9); 

Time: 

(D16) [0.384d0] 

(C17) laplace(sin(a*t)-a*t*cos(a*t),t,s); 

2 

A 2 S 1 

(D17) A ( --) 

2 2 2 22 2 2 
S+A (S+A) S+A 

(C18) time(dl7); 

Time: 

(D18) [0.233d0] 

%Inverting the Laplace transform. 

(C19) ilt((s+6)/(s A 2+4*s+12),s,t); 

-2T 2 SIN(2 SQRT(2) T) 

(D19) %E ( + COS(2 SQRT(2) T» 

SQRT(2) 

(C20) time(d4); 

Time: 

(D20) [0.933d0] 


%Transforming a differential equation into algebraic equation. 
(C21) laplace(diff(y(x),x,2)-3*diff(y(x),x)+2*y(x)=0,x,s); 


d ! 

(D21) (Y(X)) ! - 3 (S LAPLACE(Y(X), X, S) - Y(0)) 

dX ! 

!X = 0 

2 

+S LAPLACE(Y (X),X,S)+2 LAPLACE(Y(X),X,S)-Y(0) S=0 

(C22) time(d21); 

Time: 
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(D22) [0.1 5dO] 

(C23) laplace(diff(y(x),x,2)+w A 2*y(x)-b*sin(w*x)=0,x,s); 


d i 2 2 

(D23)- — (Y(X)) ! +W LAPLACE(Y (X),X,S)+S LAPLACE(Y(X),X,S) 
dX ! 

!X=0 

B W 

Y(0) S = 0 

2 2 
W +S 

(C24) time(d23); 

Time: 

(D24) [0.217d0] 


(c) Series expansion — The MACSYMA version in University of Michigan provides the 
Taylor series and power series expansion capabilities. Although the function of Fourier 
series expansion is available on the market, it is not available here. 

(C25) taylor(%e A x,[x,0,7]); 

2 3 4 5 6 7 

X X X X X X 

(D25)/T/ 1 + X + - + — + — + — + + + . . . 

2 6 24 120 720 5040 

(C26) time(d25); 

Time: 

(D26) [0.567d0] 

(C27) powerseries(%e A x,x,0); 

INF 

=== II 

\ X 

(D27) > - — 

/ 11 ! 


11=0 

(C28) time(d27); 

Time: 

(D29) [0.45d0] 
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CHAPTER IV 


APPLICATION OF SYMBOLIC AND ALGEBRAIC MANIPULATION TO 
AUTOMATIC PROBLEM FORMULATION 


4.1 Introduction 


e the most important advantages gained from symbolic and algebraic manipulation 
to automatically formulate lengthy mathematical equations without making any enont 
odem scientists and engineers are continually being challenged with mom and more 
complicated fonnulas. With the aid of symbol, c and algebmic manipulate, most problems can 
be treated eastly and correctly. Thts chapter will demonstrate six automatic formulation 
examples done by symbolic and algebraic manipulation. They are : 

1. The derivation of equations of motion in dynamics. 

2. Tensor formulation for the shell problem. 

3. The approximation to a function by Fourier series. 


4. The formulation template for the iteration method in nonlinear numerical analysis. 

5. Finite element stiffness matrix and mass matrix construction for 6-node triangular 
element in a heat transfer problem. 


6. Finite element stiffness matrix construction for 4-node isoparametrically quadrilateral 
element in a plane elasticity problem. 


Of course, the use of symbolic and algebraic mantpulators as tools to automatically 
omtulate mathematical equauons can be extended to any fields. Although the methodologies 

staitT 0n ' he Pr0b ' em ‘° * SO ' Ved ' ‘ he t>aS ' C COmma " ds “““ m programming are 
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4.2 Derivation of equation of motion by SAM 
4.2.1 Introduction 

The derivation of equations of motion by the Lagrange method for the system shown in 
Ftgure 4. 1 tnvolves finding kinetic energy T, potential energy V, and therefore the Lagrangtan 



The Lagrangian then is partially differentiated with respect to both general, zed 
coordinates and the rate of generalized coordinates. The equations of motion will be obtained 
after taking the time derivatives to appropriate terms and assembling the necessary terms. 
Mathematically, the Lagrange equations are expressed in the form of 


where L is the Lagrangian and is defined as the difference between kinetic energy T and 
potential energy V. c ' 


T - * ri)‘ * TM l( f») ! ♦ iW (4 2) 

V - - M.grfl - cos 8) - (r - r cos 8 * x sin $) * k* ( r ff)‘ (4.3) 


4.2.2 REDUCE program and solution 

1: on time; 

Time: 83 ms 
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%declaring theta and x as a function of time. 
2: depend theta,time; 

Time: 67 ms 

3: depend x.time; 

Time: 50 ms 


%Calculating the total kinetic energy of system. 

%The velocity v of body 2 need be evaluated later. 

4: te:=(l/2)*i0*(df(theta,time))**2+(l/2)*ml*r**2*(df(theta,time))**2+(I/2)*m2*v**2; 

2 2 2 2 
DF(THETA,TIME) *10 + DF(THETA,TIME) *R *M1 + M2* V 


TE := 

Time: 500 ms 


%Calculating the position vector of body 2. 

5: p:=(r*cos(theta)-x*sin(theta))*j+(r*sin(theta)+x*cos(theta))*i- 
t - C OS(T H ET A )*I*X + COS(T H ETA)*J*R + SIN(THETA)*I*R-SIN(THErA)*J*X 

x 1 1 iie* xxx $ 


%The velocity vectors is then obtained by taking the derivative of 
%the position vector with respect to time. 

6: dp:=df(p,time); 


DP := COS(THETA)*DF(THETA,TIME>*I*R-COS(THETA>* DRTHETA 
TIME)*J*X+COS(THETA)*DF(X,TIME)*I-DF{THErA TIME)* 

*sSIt^^** 7 x ^ ' DF ^ THETA,TI ^ ) * SWHETA) * jiR * DF ^ x ’ TIME) 

Time: 234 ms 


% Getting the magnitude square of velocity of body 2. 

7: v* * 2: =lcof(dp,i)* *2+lcof (dp j)* *2; 

2 2 2 2 2 2 
V :-COS(THETA) *DF(TH ETA .TIME) *R +COS(THETA) * DF(THET A .TIME) 

2 2 2 
*X +2*COS(THETA) *DF(THETA,TIME)*DF(X,TIME)*R4COS(THEFA) * 

2 2 2 2 2 
DF(X.TIME) +DF(THETA,TIME) *SIN(THETA) *R +DF(THETA .TIME) * 

2 2 

SIN(THETA^*X +2*DF(THErA,TIME)*DF(X,TIME)*SIN(THETA)*R+ 

DF(X,TIME) *SIN(THETA) 

Time: 650 ms 


%Teaching REDUCE the trigonometric rule for simplification of %formulae 
8: let (cos(theta))** 2 +(sin(theta))** 2 =l; 

Time; 150 ms 
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%Inputting the potential energy of system. 

9: ve:=-ml*g*r*(l-cos(theta))-m2*g*(r*(l-cos(theta))+x*sin(theta))+(l/2)*k*(r*theta)**2- 
VE := (2*COS(THETA)*G*R*Ml+2*COS(THETA)*G*R*M2-2*SIN(THETA) 

2 2 

*G*M2*X- 2*G*R*M1 - 2*G*R*M2 + K*R *THETA )/2 
Time: 483 ms 


%calculating the Lagrangian. 

10: la=te-ve; 

LA := -(2*COS(THETA)*G*R*Ml+2*COS(THETA)*G*R*M2-DF(THETA, 
2 2 2 2 2 
TIME) *10 - DF(THETA,TIME) *R *M1 - DF(THETA,TIME) *R 

*M2 - DF(THETA,TIME) *M2*X -2* DF(THET A ,TI ME) *DF(X ,TI ME) 

*R*M2 - DF(X,TIME) *M2 -2*SIN(THETA)*G*M2*X-2*G*R*M1- 
2 2 

2*G*R*M2+K*R *THETA )/2 
Time: 617 ms 


%Deriving the equation of motion for theta coordinate without 
%any simplification. 

1 1: el:=df(df(la,df(theta,time)),time)-df(la, theta); 

El := -(COS(THETA)*G*M2*X+DF(THETA,THETA,TIME)*DF(THETA, 

2 

TIME) *IOfDF(THETA,THETA,TIME)*DF(THETA,TIME)*R *M1 

2 

+DF(THETA, THETA, TIME)*DF(THETA,TIME)*R *M2 +DF(THETA 

2 

THET A ,T I ME) * DF(THET A .TIME) * M2* X +DF(THETA .THETA .TIME) 

2 

*DF(X,TIME)*R*M2-DF(THETA,TIME,2)*I0-DF(THETA,TIME,2)*R 

2 2 
*M1-DF(THETA,TIME,2)*R *M2-DF(THETA,TIME,2)*M2*X -2* 

DF(THETA,TIME)*DF(X,TIME)*M2*X-DF(X,TIME,2)*R*M2 + 

2 

SIN(THETA) *G*R*M1+SIN(THETA)*G*R*M2-K*R *THETA) 

Time: 783 ms 


^Deriving the equation of motion for x coordinate without simplification. 

12: e2:=df(df(la,df(x,time)),time)-df(la,x)+f ; 

2 

E2 := DF(T HET A ,TI ME, 2) * R* M2- DF(THET A ,TI ME) *M2*X-DF(THETA 
TIME) *DF(X,TIME,X)*R*M2-DF(X,TIME,X)*DF(X,TIME)*M2 + ’ 
DF(X,TIME,2)*M2 - SIN(THETA)*G*M2+F 
Time: 483 ms 


%Teaching REDUCE the 2nd derivative of theta w/t theta and time.is null, and so does x. 
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13: let df(theta,time,theta)=0,df(x,x,time)=0,df(theta, theta, time)=0,df(x,time,x)=0; 

Time: 350 ms 

%tuming off the automatic expansion switch. 

14: off exp; 

Time: 50 ms 


%Starting polynomial manipulation to simplify the equations of motion. 


%Getting the coefficient of angular acceleration term. 

15: al:=lcof(el,df(theta,time,2)); 

2 2 
A 1 := (Ml + M2)*R + 10 + M2*X 
Time: 617 ms 

%Getting the coefficient of sin(theta). 

16: a2:=lcof(el,sin(theta)); 

A2 := - (Ml + M2)*G*R 
Time: 550 ms 

%Getting the leftover after taking off the above two terms. 

17: a3:=el-al*df(theta,time,2)-a2*sin( theta); 

A3 :=-(COS(THETA)*G*M2*X-2*DF(THETA,TIME)*DF(X,TIME)*M2*X 

2 

-DF(X,TIME,2)*R*M2-K*R *THETA) 

Time: 533 ms 

%Rearranging the equation of motion for theta coordinate. 

%The results are simpler than those of in command 1 1. 

18: el:=al*df(theta,time,2)+a2*sin(theta)+a3; 

2 2 

El := ((M1+M2)*R +I0+M2*X )*DF(THETA,TIME,2MCOS(THETA)*G*M2* 

2 

X-2*DF(THETA,TIME)*DF(X,TIME)*M2*X-DF(X,TIME,2)*R*M2-K*R 
*THETA)-(M1+M2)*SIN(THETA)*G*R 
Time: 400 ms 

^Getting the coefficient of x acceleration term. 

19: cl:=lcof(e2,df(x,time,2)); 

Cl := M2 
Time: 350 ms 

%Getting the coefficient of angular acceleration term. 

20: c2:=lcof(e2,df(theta,time,2)); 

C2 := R*M2 
Time: 333 ms 

%Getting the leftover terms by removing the above two terms. 

21: c3:=e2-cl*df(x,time,2)-c2*df(theta,time,2); 

2 

C3 := - (DF(THETA.TIME) *M2*X + SIN(THETA)*G*M2 - F) 
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Time: 350 ms 


%equation of motion for x coordinate. 

22: e2:=c 1 *df(x, time, 2)+c2*df(theta, time, 2)+c3; 

2 

E2 := -((DF(THETA,TIME) *M2*X+SIN(THETA)*G*M2-F)- 
DF(THET A ,TI ME,2)* R* M2 - DF(X,TIME,2)*M2) 

Time: 317 ms 

23: bye; 

The examples shown above are the demonstration of obtaining the left hand side 
formulae of equation (4.1). The equations of motion of system can be simply done by setting 
the results of command 18 and 22 equal to zero. As the system is complex, the analytical 
derivations of equations of motion by hand will become tedious and prone to error. For the 
cases of complex systems. The application of symbolic and algebraic manipulation will be 
more significant. 

4.3 Automatic tensor formulation for shell problem 
4.3.1 Preliminary formulation 


Tensors are convenient mathematical entities for concisely describing physical 
situations, that are independent of coordinate transformations. Although the benefits gained by 
employing tensor notation are quite significant in the related fields, the expressions of tensor 
formula often rise to errors. Fortunately, this difficulty can be avoided by the use of symbolic 
and algebraic manipulation. This advantage will be demonstrated by formulating the thin shell 
problem in tensor form and using SAM to expand the resulting tensor equations. For the sake 
of convenience, the Latin indices will be referred for the range 1, 2, 3 and the Greek indices are 
in the range of 1, 2 in the following paragraph without special stress. 


On the formulation of the thin shell problem, the only necessary inputs are three 
i 1 2 

parametric equations f (u ,u ) of the shell middle surface. Based on parametric equations, the 
covariant metric tensor a . and its determinant are calculated 


a -/7‘ .ML.JL 

“aft J ,f> du a du f> 
a - det (a^ ) 


(4.4) 

(4.5) 


The contravariant metric tensor is therefore given by 
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X* 



a* - cofactor ( a ^ )la (4 6) 

The first and second kind of Christoffel symbols for the surface are defined as follows 
respectively : 

J - 7<a* , a -a^ y ) (4.7) 

{/9°)'}** a ‘' ( ‘V. +<, <, .»-%<) (48) 

The covariant differentiations of a covariant vector (the 1st order tensor) and the 2nd order 
tensor are therefore calculated 

V«- D « A > “ Va“{/ a } A r (4.9) 

A »-- D ° A » - {/„}.«* - { r ( 410 ) 

When the above formulations are done, the curvature tensor and its determinant are then 
computed by the formulae of 
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(4.11) 

(4.12) 


d -det(rf^) 


Where the X 1 are the unit normal vectors of the shell middle surface [see Figure 4.2] and are 
equal to 


X' 



(4.13) 


Where is the generalized Kronecker delta. 

The strain tensor is defined as half of the difference between the deformed metric tensor 
and undeformed metric tensor. 


“ 2 ( a a0 a c# ) 


(4.14) 


Where the asterisk superscript is denoted as deformed state. The a* a p is defined as 

^ + P* + P ^ + P 1 aP, 1 ^q a q it (4.15) 


The generalized two-dimensional displacement gradient p, and rotation q are represented as 

P* -D a v p - d^w (4.1 6) 

4° -d^v* + w a (4.17) 


Where the v Q are the in-plane displacements and w is the out of plane displacement of 
the middle surface of the shell. Similarly the bending tensor is equal to the difference between 
the deformed curvature tensor and undeformed curvature tensor. 


K 


<*> 




(4.18) 


The deformed curvature tensor d a p is expressed as 

a . i/2 r 


- <F>' 10 + *D t q,*d\p„) 




Where =det (p a p) 

The constitutive equations for isotropically elastic and thin shell are 
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(4.19) 


N <* - — -£2L_ [(1 _ y)E < * + w*E r r ] 

1 - v r 

M *~ i2 ( f- 1(1 ~ v) ** 1 < 4 - 2 °) 

afi 

Where the Young's modules E and strain tensor E shouldn't be confused here.The 
final equilibrium equations are then 

DaN* + 2dy D a MT + D a dy + F* - 0 ( 4 . 21 ) 

DaDpM 0 * -P - 0 (4.22) 

Where the and P are the external forces applied in the in-plane and out-plane directions, 
respectively. 


4.3.2 REDUCE program and resultant expressions 


The REDUCE program shown below is based on the above formulation methodology. 
The explanations of the program are also included to facilitate an understanding where it is 
necessary. The resultant expressions are too huge to be included here and are available in 
reference [17]. 

ARRAY X(3),C1(2,2,2),C2(2,2,2); 

OPERATOR U,V,W,F; 

MATRIX A(2,2),CONTRA(2,2),D(2,2) ; 

% INPUTTING SURFACE PARAMETRIC FUNCTIONS 

X(1):=U(1); 

X(2):=U(2); 

X(3):=CONSTANT ; 

%=========== = === = == =========================== _ ======= 

% CALCULATING COVARIANT METRIC TENSOR & ITS DETERMINANT 

FOR M:=l:2 DO FOR N:=l:2 DO 

A(M,N):=FOR I:=l:3 SUM DF(X(I),U(M))*DF(X(I),U(N)); 

DETA:=DET(A); 

DEPEND W,U(1),U(2); 

FOR I:=l:2 DO DEPEND V(I),U(1),U(2); 

%CALCULATING CONTRA VARIANT METRIC TENSOR COMPONENT* 
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FOR L:=l:2 DO FOR M:=l:2 DO 
IF L=1 and M=1 THEN CONTRA(L,M):=A(2,2)/DETA 
ELSE IF L NEQ M THEN CONTRA(L,M):=-A(M,L)/DETA 
ELSE CONTRA(L,M):=A( 1,1)/DETA; 

%CALCULATING THE 1ST CHRI ST OLFFEL SYMBOL 


FOR L:=l:2 DO FOR M:=l:2 DO FOR N:=l:2 DO 

C1(L,M,N):=(1/2)*(DF(A(L,N),U(M))+DF(A(M,N),U(L))-DF(A(L,M),U(N))); 
%CALCULATING THE 2ND CHRISTOLFFEL SYMBOL (SEE Iq746) 


FOR L:=l:2 DO FOR M:=l:2 DO FOR N:=l:2 DO 
C2(L,M,N):=FOR I:=l:2 SUM A(L,I)*C1(M,N,I); 


%SUBROUTI NE FOR CALCULATING THE COVARIANT DERIVATIVE 
%FOR THE 1ST ORDER TENSOR 

PROCEDURE COVD(L,VAR(M)); 

DF(VAR(M),U(L))-(FOR I:=l:2 SUM C2(I,L,M)*VAR(I)); 

^SUBROUTINE FOR CALCULATINE THE COVARFAt^lS^ _ VATrvi _ 
%FOR THE 2ND ORDER TENSOR 

PROCEDURE COVD2(L,FUN(M,N)); 

DF(FUN(M,N),U(L))-(FOR I:=l:2 SUM C2(I,L,M)*FUN(I,N)) 

-(FOR I : = 1 : 2 SUM C2(I,L,N)*FUN(M,I)); 

MATRIX E(2,2),K(2,2),P(2,2),FF(2 > 3),XX(1,3); 

%CALCULATING RELATIVE ALTERNATE TENSOR 
%ALSO CALLED AS GENERALIZED KRONECKER DELTA 

PROCEDURE EE(I,J,K); 

IF I=J OR J=K OR K=I THEN 0 
ELSE IF 

(1=1 AND J=2 AND K=3) 

OR (1=3 AND J=1 AND K=2) 

OR (1=2 AND J=3 AND K=l) 

THEN 1 
ELSE -1; 


%CALCULATING THE CURVATURE TENSOR & ITS DETERMINANT 
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FOR I:=l:2 DO FOR J:=l:3 DO 
FF(I,J):=DF(X(J),U(I)); 

FOR I:=l:3 DO 

XX( 1 ,1):=( l/SQRT(DETA))*(FOR J:=l:3 SUM 

(FOR K:=l:3 SUM EE(U,K)*FF(1,J)*FF(2,K))); 

FOR I:=l:2 DO FOR J:=l:2 DO 
D(I,J):=FOR 1 1:=1:3 SUM XX(1,I1)*DF(FF(I,I1),U(J)); 
DETD:=DET(D); 

^CALCULATING THE GENERIZED DISPLACEMENT tInSOR ~” 


FOR L:=l:2 DO FOR M:=l:2 DO 

P(L,M):=DF(V(M),U(L))-(FOR I:=l:2 SUM C2(I,L,M)*V(I))-D(L,M)*W; 

%CALCULATING THE ROTATION TENSOR (SEE EQ. 4.17) 

ARRAY Q(2); 

FOR M:=l:2 DO 

Q(M):=DF(W,U(M))+(FOR I:=l:2 SUM D(M,I)*(FOR J:=l:2 SUM CONTRA(I,J)*V(J))); 
%DERI VING THE STRAIN TENSOR 


FOR I:=l:2 DO FOR J:=l:2 DO 
E(I,J):=(1/2)*(P(I,J)+P(J,I) 

+(FOR Il:=l:2 SUM (FOR Jl:=l:2 SUM CONTRA(Il,Jl)*P(I,Jl))* P(J,I1))+Q(I)*Q(J)) 
OFF PERIOD; 

ON FORT ; 

OFF PERIOD; 

FOR I:=l:2 DO FOR J:=l:2 DO 
WRITE - E( H ,I,", H ,J,")= M > E(I,J); 


%CALCULATING THE ABSOLUTE 2-D ALTERNATE TENSOR 
%============== ============================= __ 


PROCEDURE EPS(L,M); 

IF L=1 AND M=2 THEN SQRT(DETA) 

ELSE IF L=2 AND M=1 THEN -SQRT(DETA) 

ELSE 0; 

%CALCULATE CONTRA VARIANT ABSOLUTE 2-D ALTERNATE TENSOR 


PROCEDURE CEPS(L.M); 

IF L=1 AND M=2 THEN 1/SQRT(DETA) 

ELSE IF L=2 AND M=1 THEN -1/SQRT(DETA) 
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ELSE 0; 


%============================ 

%CALCULATING THE BENDING TENSOR 
%============================== 


FOR I:= 1:2 DO FOR J:=l:2 DO 

«K(I,J):=SQRT(DETA/AD)*((l+(FOR L:=l:2 SUM CONTRA(L,L)*P(L,L)) 
+DET(P)/DETA)*(D(I,J)+DF(Q(I),U(J)) 

-(FOR II := 1:2 SUM C2(I1,J,I)*Q(I1)) 

+(FOR L:=l:2 SUM (FOR M:=l:2 SUM CONTRA(L,M)*D(J,M))*P(I,L))) 
-(FOR L:=l:2 SUM ((FOR M:=l:2 SUM CONTRA(L,M)*Q(M)) 

+(FOR N:=l:2 SUM CEPS(L,N)*(FOR S:=l:2 SUM 

(FOR R:=l:2 SUM CEPS(R,S)*Q(R))*P(S,N))))*(DF(P(I,L),U(J)) 

-(FOR 12:— 1:2 SUM C2(I2,J,I)*P(I2,L)) 

-(FOR I2:= 1:2 SUM C2(I2,J,L) 

*P(I,I2))-D(J,L)*Q(I))))-D(I,J); 

WRITE " K( H ,i, H ,",j,")=",K(I,J)»; 

DEPEND AD,U(1),U(2); 

RAD:=DET(A+2*E)$ 

WRITE " AD=",RAD; 


%=================================== 

%DERI VING THE CONSTITUTIVE EQUATIONS 
%================================ 

MATRIX N(2,2),M(2,2); 

FOR I:=l:2 DO FOR J:=l:2 DO 


%- - — 

%EFFECTIVE MEMBRANE STRESS TENSOR 
% 

«N(U):=(YE/( 1-V V**2))*(( 1- V V)*(FOR L:=l :2 SUM CONTRA(L,I)* 
(FOR 1 1:=1:2 SUM CONTRA(1 1,J)*E(L ( I l)))+VV*CONTRA(U)* 
(FOR L:=l:2 SUM CONTRA(L,L)*E(L ) L))) ; 

% 

%EFFECTIVE MOMENT TENSOR 

% 

M(U):=(YE*H**3/(12*(1-VV**2)))*((1-VV) 

*(FOR L:=l:2 SUM CONTRA(L,I)* 

(FOR Il:=l:2 SUM CONTRA(IU)*K(L,Il))) 
+VV*CONTRA(I,J)*(FOR L:=l:2 
SUM A(L,L)*K(L,L)))»; 


%=========================—======================== 

% WRITING-OUT EXPRESSIONS OF STRESS AND MOMENT TENSOR 
%=================================================== 


FOR I:=l:2 DO FOR J:=l:2 DO 
WRITE " N('\I,V,J,")=",N(I,J); 

FOR I:=l:2 DO FOR J:=l:2 DO 



WRITE " 


%DERIVING THE EQUILIBRIUM EQUATIONS 

%======================~==== ====== 


ARRAY M1(2),N1(2); 

FOR I:=I:2 DO 
Ml(I):=FOR J:=l:2 SUM 

(DF(M(I,J),U(J))+(FOR 1 1:=1:2 SUM C2(I,I1,J)*M(I1,J)) 
+(FOR 1 1:=1:2 SUM C2(J,I1,J)*M(I,I1))); 

MM:=FOR I:=l:2 SUM 

(DF(Ml(I),U(I))+(FOR 1 1:=1:2 SUM C2(I,I1,I)*M1(I1))); 
BEQ:=MM-(FOR Il:=l:2 SUM 
(FOR I2:= 1:2 SUM D(Il,I2)*(FOR I3:=l:2 SUM 
CONTRA(I2,I3)*(FOR I4:=l:2 SUM D(I4,I3)*M(I1,I4))))) 
-(FOR 1 1:= 1:2 SUM 

(FOR I2:=l:2 SUM D(I1,I2)*N(I1,I2)))-PP; 

FOR I:=l:2 DO 
Nl(I):=FOR J:=l:2 SUM 

DF(N(I,J),U(J))+(FOR 1 1 := 1:2 SUM C2(I,I1,J)*N(I1,J)) 
+(FOR Il:=l:2 SUM C2(J,I1,J)*N(I,I1)) 

+2*(FOR Il:=l:2 SUM CONTRA(I,1 1)* 

(FOR I2:=l:2 SUM D(I2,U)*(DF(M(I2,J),U(J)) 

+(FOR I3:= 1:2 SUM C2(I2,i3,J)*M(I3,J)) 

+(FOR 13: =1:2 SUM C2(J,I3,J)*M(I2,I3)))))- 
(FOR Il:=l:2 SUM M(I1,J) 

*(FOR I2:=l:2 SUM D(I1,I2)*(DF(C0NTRA(I,I2),U(J)) 
+(FOR I3:=l:2 SUM C2(I,I3,J)*CONTRA(I3,I2)) 

-KFOR I3:=l:2 SUM C2(I2,I3 ) J)*CONTRA(I,I3)))) 

+(FOR I2:=l:2 SUM CONTRA(Il,I2)*(DF(D(IU2),U(J))- 
(FOR I3:=l:2 SUM C2(I3,I1,J)*D(I3,I2))- 
(FOR I3:=I:2 SUM C2(I3,I2,J)*D(1 1,I3)))))+F(I); 

WRITE BEQ," =0"; 

FOR I:=l:2 DO WRITE N1(I)," =0"; 

BYE; 


4.3.3 Remarks 


The REDUCE program and solution shown above are just for the case of plates. For 
the geometry other than plates, the program is easily modified by changing the input parametric 
equations in the beginning of program. As the results show, the formulae of strain, bending, 
stress, moment tensors and equilibrium equations are functions of three displacements and their 
derivatives. Therefore we may assume three displacement fields as polynomial (or power 
series, Fourier series etc.), neglect the unnecessary terms and solve the problem. Since it is out 
of the scope of this report, it will be left to the interested researchers. 
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4.4 Approximation of a function by Fourier series 
4.4.1 Introduction 

Mathematically, two functions which belong to different spaces are apparently different 
in properties. It is theoretically impossible to replace one in terms of the other. However, to 
approximate one by the other is feasible and is actually adopted by many engineers in various 
fields. The purpose of approximation depends on the problem. It may be for avoiding the 
mathematical difficulties or for simplifying the function so that it can be solved easily. Some 
popular methods to approximate a function include Taylor series approximation and Fourier 
series expansion etc. Since they are the approximation methods, the accuracy is up to the 
number of terms included. Of course, as more terms are involved in the formulation, more 
complicated expression will turn out and possibility of making errors is increased. Based on 
these reasons, even though a function can be approximated theoretically, a high accuracy in 
evaluation of such approximations is difficult to achieve. By the symbolic and algebraic 
manipulation, the approximation can be formulated without errors and the desired accuracy can 
be attained. Here, an example of Fourier series approximation to the hat function will be shown 
to demonstrate the advantages of application of symbolic and algebraic manipulation. 

The hat function is a fundamentals of shape function in finite element method. It is 
defined as : 

f/ {x ) - 1 + x , for - 1 s x < 0 

[/ (x ) - 1 - x , for 0 s x < - 1 

Mathematically, it's a piecewise C 1 continuous function which certainly is different 
from the sinusoidal function which is C 00 continuous function within the domain they span. 
The approximation to the hat function by Fourier series is expressed as 

/(* ) “ — + 2, anCOS (~T~ > + ) ( 4 . 23 ) 

~ n -1 

where the Fourier coefficients are evaluated as follows 
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(4.24) 


a o - 7 -/ / (* )dx 
a n - i - //^ ) cos (—^—)dx (4.25) 

b n ~ J^f j( x )sin(^-)dx (4.26) 


4.4.2 REDUCE program for generating Fourier series 

The REDUCE program to generate the Fourier series and output a fortran subroutine 
for hat function is shown in the follow assuming L=1.0. This program can be used to generate 
Fourier series for an arbitrary function by simply changing the input function. 


% Inputting the given function and informations. 

% M : number of piecewise bounded interval. 

% f(M): the function in the Mth interval. 

% c(n) : the upper and lower limit of finite integration. 
% 1 : half length of interval. 

% k : number of Fourier series terms needed. 


M:=2; 

K'=50’ 

ARRAY F(M),C(M+1); 
f(l):=l+x; %c(l) =< x =< c(2) 
f(2):=l-x; %c(2) =< x < c(3) 
c( 1 ): =- 1 ;c( 2) : =0;c(3): = 1 ; 
l:=(c(m+l)-c(l))/2; 

% Obtaining the Fourier series coefficients 


aO:=for i:=l:m sum 

(sub(x=c(i+l),int(f(i),x))-sub(x=c(i),int(f(i),x)))/l; 
an:=for i:=l:m sum 

(sub(x=c(i+l),int(f(i)*cos(n*pi*x/l),x)) 
-subfx^O.intfffO^cosf^pi^xyO.x)))/!; 
bn:=for i:=l:m sum 

(sub(x=c(i+l),int(f(i)*sin(n*pi*x/l),x)) 
-sub(x=c(i),int(f(i)*sin(n*pi*x/l),x)))/l; 
on rat; 
on div; 

%========= =========== 

% Generating the Fourier series. 
%============ ======== 
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fs:=aO/2+for i:=l:k sum 

sub(n=i,an)*cos(i*pi*x/l)+sub(n=i,bn)*sin(i*pi*x/l); 
off echo; 
on fort; 
cardno!*:=10; 


%=========— ======================= 

% Outputting the fortran subroutine of Fourier series. 
%================================== 

out "fourier.ftn"; 

write " subroutine fourier(x,fs)"; 
write " implicit real*8(a-h,o-z)"; 
write" pi=3. 141592654"; 
fs:=fs; 

write" return"; 
write" end"; 
shut "fourier.ftn"; 
bye; 


4.4.3 Resultant fortran subroutine from REDUCE 


The following results are produced automatically from the above REDUCE program for 
fifty terms Fourier series case. This subroutine can be directly input into fortran main program 
without any troubles. 


subroutine fourier(x,fs) 
implicit real*8(a-h,o-z) 
pi=3. 141592654 

ANSl=4./625.*COS(25.*PI*X)*PI**(-2)+4./529.*COS(23.* 

PI*X)*PI**(-2)+4./441.*COS(21.*PI*X)*PI**(-2)+4./ 

. 361.*COS(19.*PI*X)*PI**(-2)+4./289.*COS(17.*PI*X)*PI 
. **(-2)+4./225.*COS(15.*PI*X)*PI**(-2)+4./169.*COS( 

. 13.*PI*X)*PI**(-2)+4./121.*COS(l l.*PI*X)*PI**(-2)+4./ 

. 81.*COS(9.*PI*X)*PI**(-2)+4./49.*COS(7.*PI*X)*PI**(-2 
. )+4./25.*COS(5.*PI*X)*PI**(-2)+4./9.*COS(3.*PI*X)*PI 
. **(-2)+l./2. 

FS=4.*COS(PI*X)*PI**(-2)+4./2401.*COS(49.*PI*X)*PI**( 
. -2)+4./2209.*COS(47.*PI*X)*PI**(-2)+4./2025.*COS(45. 

*PI*X)*PI**(-2)+4./1849.*COS(43.*PI*X)*PI**(-2)+4./ 

. 1681.*COS(41.*PI*X)*PI**(-2)+4./1521.*COS(39.*PI*X)* 

. PI**(-2)+4./ 1369. *COS(37.*Pl*X)*PI**(-2)+4./ 1225.* 

. COS(35.*PI*X)*PI**(-2)+4./1089.*COS(33.*PI*X)*PI**( 

. -2)+4./961.*COS(31.*PI*X)*PI**(-2)+4./841.*COS(29.* 

. PI*X)*PI**(-2)+4./729.*COS(27.*PI*X)*PI**(-2)+ANSl 

return 

end 
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Four curves are plotted in the same figure for three-term, five-term , seven-term and 
fifty-term cases. As the figure shows, the hat function can be actually simulated by Fourier 
series. When fifty terms are used, the difference between hat function and its Fourier series is 
just invisible although they are in different spaces. This is one of the advantages of the 
application of symbolic and algebraic manipulation. 



Figure 4.3 : Convergence of Fourier series approximation 


4.5 Template for nonlinear numerical analysis 
4.5.1 Introduction 

In nonlinear numerical analysis, it is necessary to evaluate the Jacobian matrix and to 
solve the system of equations at each iteration. Symbolically, the Newton-Raphson iteration 
can be expressed as 


(4.27) 

Since the evaluation of the inverse of the Jacobian matrix is quite time-consuming, 
equation (4.27) is traditionally changed into the following form and then is solved as a linear 
system of equations at each iteration stage. 

[J]{AX} < ‘ ) --{F} (4.28) 

Although the avoidance in evaluation of the inverse of the Jacobian does expedite the 
execution, it still takes time to solve a system of equations, especially for the case of a large 
number of unknowns. With the help of symbolic and algebraic manipulation, the efficiency can 
be improved further. 
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Instead of transforming equation (4.27) into equation (4.28), it is rearranged into 

{X} tk '-{xf*-{LX} tki --[jr l {F} (4.29) 

The application of symbolic and algebraic manipulation to this problem is to make a 
template form of the right, hand, side of equation (4.29). This includes the symbolic evaluation 
of the Jacobian matrix, its inverse, and the multiplication of the matrix by load vector. The 
results are then converted into fortran code for numerical analysis. The iterations are done by 
simply substituting the current solutions into the template formulae to get the residuals. This 
substitution is much faster than solving the system of equations. This is a new way to improve 
the efficiency of program execution. 

4.5.2 Preliminary formulation 

The above idea is demonstrated by solving the Fermat- Weber location problem for the 
case of a two-dimensional Euclidean space. Given n fixed points .(Xj, yp, i=l,..,n, in the 
plane, the object of the Fermat- Weber problem is to find the best location to minimize the total 
length of the distances among the optimal location and each point. Mathematically, it is 

Minimize -x t ) 2 + (y -y ) 2 (4.30) 


or rewritten in the following form 

n 

Minimize 

(Si 1 2 

where (x,y) is the coordinate of the optimal location to be found, and 



(4.31) 


(4.32) 

(4.33) 

(4.34) 


Intuitively, the minimization will be accomplished by differentiating equation (4.31) 
with respect to Y and setting the results to zero. But since the objective function in equation 
(4.31) is not differentiable at the exact (x-,y-) points, a smoothing parameter e is introduced 
into the object function to avoid the difficulty. 
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+ e 


(4.35) 


®,(Y ) 





The perturbed ^bjective function now becomes differentiable everywhere and is strictly 
convex. Then -%r - 0 will result in a fixed point iteration 


-0,1,2 


(4.36) 


Where the weight wW is defined as 


w 


<*) 




2 


(4.37) 


4.5.3 REDUCE program 


Consider the specific problem which is to find the optimal location among three points 
(0,1), (0,-1) and (x,0) where x varies in the range [0, °o). The REDUCE program is based on 
equation (4.36). 


% 

% waa : W*A*A in equation (4.36) 

% wac : W*A*C in equation (4.36) 

% 

matrix waa(2,2),wac(2,l); 

waa( 1 , 1 ) : = 1 /sqrt(r 1 )+ 1 /sqrt( r2)+ 1 /sqrt( r3) ; 

waa(2,2):=waa( 1,1); 

wac(l,l):=xl/sqrt(rl)+x2/sqrt(r2)+x3/sqrt(r3); 

wac(2,l):=yl/sqrt(rl)+y2ysqrt(r2)+y3/sqrt(r3); 

wac:=( 1 /waa)* wac; 

on fort; 

off echo; 

off period; 

cardno!*:=10; 

out "cssa.ftn"; 

write " subroutine cssa(xl > yl,x2,y2,x3,y3,pu,x0,y0)"; 

write" implicit real *8(a-h,o-z)"; 

write " rl=(x-xl)**2+(y-yl)**2+pu**2"; 

write " r2=(x-x2)**2+(y-y2)**2+pu**2"; 

write " r3=(x-x3)**2+(y-y3)**2+pu**2"; 

write" x=",wac(l,l); 

write " y=",wac(2,l); 

write " return"; 

write " end"; 

shut "cssa.ftn"; 

bye; 
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4.5.4 Fortran subroutine from REDUCE 


The fortran subroutine produced here is applicable when the number of points is three 
The coord, nates of three points are arbitrary. The perturbation parameter is an arbitranly small 
number except zero. The resultant optimal locations are plotted in Figure 4.4. with respect to 
variable x3. The execution time for this problem on Apollo workstation Domain 4000 is too 

smaU to be measured. The difference in execution time will be more significant when the size 
oi the problem is increased. 

subroutine cssa(xl,yl,x2,y2,x3,y3,pu,x,y) 
implicit real*8(a-h,o-z) 
r l=(x-xl)**2+(y-y l)** 2 +pu **2 
r 2 =(x-x 2 )** 2 +(y-y 2 )** 2 +pu **2 
r3=(x-x3)**2+(y-y3)**2+pu**2 

J^SQRT(R2)*SQRT(Rl)*R3*Xl+SQRT(R2)*SQRTfRn*R3*X?-i- 

} * S Q RT ( R 1 ) * 1 +SQRT (R3) * SQRTYR 1 )* R2^3+SORT 

^ * ^QRT(R2)*R1*X2+SQRT(R3)*SQRT(R2)*R1*X3+R1*R^*y^ 

. +Rl*R3*X2+R2*R3*Xl)/(2*SQRT(R2)*SORT(Rn*R3+?*<5nRTf' 

return 
end 



Figure 4.4 : Behavior of optimal location vs. variable x3 
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4.6 Triangular stiffness and mass matrix construction 
4.6.1 Preliminary formulation 

It is well known that the finite element method has become a powerful tool in solving 
general engineering problems. Usually the finite element method is applied in three steps: 

1 . Domain discretization (pre-processing). 

2. Constructing the stiffness, mass matrices and load vectors etc., then prescribing the 
boundary conditions and solving them. 

3. Results treatment (post-processing). 

While computing stiffness and mass matrices, partial differentiation, matrix 
multiplication, matrix inversion and integration are required. For a higher order interpolation, 
the formulation is always very tedious and prone to introduce errors. With the aid of symbolic 
and algebraic manipulation, all these troubles can be alleviated. The example shown in this 
section will demonstrate the application of symbolic and algebraic manipulation to the automatic 

construction of stiffness and mass matrix of six-node triangular element for a heat transfer 
problem. 


The stiffness K and mass M matrices are defined as follows assuming unit thickness. 
* -JB' * D * BdA ( 4 . 3 8 

M - Jn t *N * p * Cp dA (439 


where the shape functions are 


N s -4*s t *s 3 , N 6 -4*s 3 * Sl , 

" 5 * ~ + * 4 ). - i.(Af 4 + N s ), N 3 -s 3 - L(N s + N 6 ) (4.40) 

The B matrix is evaluated by the formulae of 
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dN l 

dN^ 

dN x 

*T 

dS 2 

*T 

<W 2 

<?a 2 

<w 2 

ds i 

35 2 


<w 3 

dN 3 



ds z 

*3 

dN 4 

<w 4 

<W 4 

ds x 

* 2 

dS 3 


dN s 


ds x 

ds i 

*3 

dN A 

o 

dN . 

o 

<?AT 

o 

*1 


*3 


The material property matrix D is assumed to be symmetric 


d n d 


The bj, c- here are expressed by global nodal coordinates. 


b i“ 

C 1 ” ~J~( X 3 ~ X 


C 1 m ~(x l ~ x 3 )' c 3 “ ~( X 2 ~ x j) 


The Jacobian J is equal to 


(4.41) 


(4.42) 


(4.43) 


J m (* l -x 3 )(y l -y 3 )-(y l - y 3 )( x 2 - * 3 ) - 2* arm (4.44) 

4.6.2 REDUCE program for stiffness matrix 

MATRIX NS(6,3),NS1(2,6),NS2(2,6),NS3(2,6),BB(2,6); 

ARRAY N(6),B(3),C(3); 

% 

%Inputting the shape functions 

% 

N(4):=4*S1*S2; 

N(5):=4*S2*S3; 

N(6):=4*S3*S1; 

N(l):=Sl-( 1/2)*(N(6)+N(4)); 

N(2):=S2-(1/2)*(N(4)+N(5)); 

N(3):=S3-( 1/2)*(N(5)+N(6)); 
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% 

% Calculate the Jacobian, area and b(i),c(i) defined in Eq.(4.43)& (4.44) 


JAC:=(X1-X3)*(Y2-Y3)-(X2-X3)*(Y1-Y3); 

AREA:=JAC/2; 

B( 1):=(Y2-Y3)/JAC; 

B(2):=(Y3-Y1)/JAC; 

C(1):=(X3-X2)/JAC; 

C(2):=(Xl-X3)/JAC; 

B(3):=-B(l)-B(2); 

C(3):=-C(l)-C(2); 

% ----- 


%Calculating B matrix defined in Eq. (4.41) 

%~ - 


FOR M:=l:2 DO « 
FOR I:=l:6 DO « 


NS(I 1 ):= DF ( N ( I ), SI ) ;NS ( 1 ) 2) : =DF( N ( I ), S2 ); N S(I3) : =DF(N(1),S3)- 

IF M=1 THEN BB(M,I):=FOR J:=l:3 SUM (B(J)*NS(U)) ELSE 
BB(M,I):=FOR J:=l:3 SUM (C(J)*NS(I,J))»»- 
FOR I:= 1 :2 DO «FOR J:= 1 :6 DO 


«NS1(I,J):=SUB(S1=2/3,S2=1/6,S3=1/6,BB(I J)V 
NS2(I,J):=SUB(S1=1/6,S2=2/3,S3=1/6,BB(I,J)); 
NS3(I,J):=SUB(S1=1/6,S2=1/6,S3=2/3,BB(I,J))»»; 
% 


%Given the symmetric material matrix. 

% 

MATRIX D(2,2); 

D:=MAT((K1 1,K12),(K12,K22)); 

MATRIX LO(6,6),NU(6,6),CC(6,6),SKE(6,6); 


%Obtaining the stiffness matrix. 
% - 


SKE:=(AREA/3)*(TP(NS1)*D*NS1+TP(NS2)*D*NS2+TP(NS3)*D*NS3)$ 

% Making the appropriate substitution to simplify the final expression. 

COB:=-(X3* Y2-X3* Y 1 -X2*Y3+X2*Y 1+X 1 *Y3-X1*Y2)' 

FOR I:=l:6 DO FOR J:=l:6 IX) 

^.j) J =NU(U)/(TCa 


%Outputting the resultant fortran subroutine 

% 

ON FORT; 

OFF ECHO; 

OFF PERIOD; 

OUT "sym"; 

subroutine Stiff(xl,yl,x2,y2 > x3 > y3 > dll,dl2,d22,ske) , '; 
WRITE" implicit real^^a-h.o-z)"; 

WRITE" dimension ske(6, 6)"; 

WRITE " DJ=-(X3*Y2-X3*Y1-X2*Y3+X2*Y1+X1*Y3-X1*Y2)" - 

FOR I:=l:6 DO FOR J:=l:6 IX) ’ 

IF J>=I THEN WRITE "SKE(",I,",",J,'')=",SKE(I J) 

ELSE WRITE "SKE(",I,",",J,*)=SKE(",J," > ",I/ ; )"; 
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WRITE " return"; 
WRITE" end"; 
SHUT "sym"; 

OFF FORT; 

BYE; 


4.6.3 Resultant stiffness matrix made by REDUCE 

subroutine stiff(x 1 ,y 1 ,x2,y2,x3 ,y3,d 1 1 ,d 1 2,d22,ske) 
implicit real*8(a-h,o-z) 
dimension ske(6,6) 

DJ=-(X3*Y2-X3*Y 1-X2*Y3+X2*Y 1+X1*Y3-X1*Y2) 
SKE(UMD11*Y2**2-2*D11*Y2*Y3+D11*Y3**2-2*D12*X2*Y2 
- +2*D 1 2*X2*Y3+2*D1 2*X3 * Y2-2*D 1 2*X3 * Y3+D22*X2**2-2* D22 
. *X2*X3+D22*X3**2)/(2*DJ) 

SKE(1,2)=(D11*Y1*Y2-D11*Y1*Y3-D11*Y2*Y3+D11*Y3**2-D12 
. *X1*Y2+D12*X1*Y3-D12*X2*Y1+D12*X2*Y3+D12*X3*Y1+D12* 

. X3*Y 2-2*D 12*X3* Y3+D22*X 1 *X2-D22*X 1 *X3-D22*X2*X3+D22* 
. X3**2)/(6*DJ) 

SKE(13)=-(D1 i*yi*Y2-Dl 1*Y1*Y3-D1 1*Y2**2+D1 1*Y2*Y3- 
. D12*X1*Y2+D12*X1*Y3-D12*X2*Y 1+2*D12*X2*Y2-D12*X2*Y3+ 
. D12*X3*Y 1-D12*X3*Y2+D22*X1*X2-D22*X1*X3-D22*X2**2+ 

. D22*X2*X3)/(6*DJ) 

SKE(1,4)=-(2*(D1 1*Y 1*Y2-D1 1*Y1*Y3-D1 1*Y2*Y3+D11*Y3**2 
. -D 12*X 1 * Y 2+D 12*X 1 * Y3-D 12*X2* Y 1+D 12*X2* Y3+D12*X3* Y 1+ 

. D12*X3*Y2-2*D12*X3*Y3+D22*X1*X2-D22*X1*X3-D22*X2*X3+ 

. D22*X3 * *2))/(3 * DJ) 

SKE(1,5)=0 

SKE(1,6)=(2*(D1 1*Y 1*Y2-D1 1*Y 1*Y3-D11*Y2**2+D1 1*Y2*Y3- 
.D12*X1*Y2+D12*X1*Y3-D12*X2*Y1+2*D12*X2*Y2-D12*X2*Y3+ 
. D12*X3*Y 1-D12*X3*Y2+D22*X1*X2-D22*X1*X3-D22*X2**2+ 

. D22*X2*X3))/(3*DJ) 

SKE(2, 1 )=SKE( 1 ,2) 

SKE(2,2)=(D1 1*Y 1**2-2*D1 1*Y1*Y3+D1 1*Y3**2-2*D12*X1*Y 1 
. +2*D12*X1*Y3+2*D12*X3*Y 1-2*D12*X3*Y3+D22*X1**2-2*D22 
. *X1*X3+D22*X3**2)/(2*DJ) 

SKE(23)=(D1 1*Y 1**2-D11*Y1*Y2-D1 1*Y 1*Y3+D1 l*Y2*Y3-2* 

. D12*X1*Y I+D12*X1*Y2+D12*X1*Y3+D12*X2*Y1-D12*X2*Y3+ 

. D12*X3*Y 1-D12*X3*Y2+D22*X1**2-D22*X1*X2-D22*X1*X3+ 

. D22*X2*X3)/(6*DJ) 

SKE(2,4)=-(2*(D1 1*Y 1*Y2-D1 1*Y1*Y3-D1 1*Y2*Y3+D1 1*Y3**2 
. -D12*X1*Y2+D12*X1*Y3-D12*X2*Y 1+D12*X2*Y3+D12*X3*Y 1 + 
.D12*X3*Y2-2*D12*X3*Y3+D22*X1*X2-D22*X1*X3-D22*X2*X3+ 

. D22* X3 * *2) )/(3 * D J) 

SKE(2,5)=-(2*(D 1 1 *Y 1 **2-D 1 1* Y 1 *Y2-D1 1 * Y I *Y3+D1 1 *Y2* Y3 
. -2*D12*X1*Y1+D12*X1*Y2+D12*X1*Y3+D12*X2*Y1-D12*X2*Y3 
. +D12*X3*Y 1-D12*X3*Y2+D22*X 1 **2-D22*X I *X2-D22*X 1 *X3+ 

. D22*X2*X3))/(3*DJ) 

SKE(2,6)=0 

SKE(3,1)=SKE(1,3) 

SKE(3,2)=SKE(2,3) 

SKE(3,3)=(D11*Y1**2-2*D11*Y1*Y2+D11*Y2**2-2*D12*X1*Y1 
. +2*D12*X1*Y2+2*D12*X2*Y 1-2*D12*X2*Y2+D22*X1**2-2*D22 
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. *X 1 *X2+D22*X2**2)/(2*DJ) 

SKE(3,4)=0 

SKE(3 ) 5)=-(2*(D11*Y1**2-D11*Y1*Y2-D11*Y1*Y3+D11*Y2*Y3 
. -2*D12*X1*Y1+D12*X1*Y2+D12*X1*Y3+D12*X2*Y1-D12*X2*Y3 
. +D12*X3*Y1-D12*X3*Y2+D22*X1**2-D22*X1*X2-D22*X1*X3+ 

. D22*X2*X3))/(3*DJ) 

SKE(3,6)=(2*(D1 1*Y 1*Y2-D1 1*Y 1*Y3-D1 1*Y2**2+D1 1*Y2*Y3- 
. D12*X1*Y2+D12*X1*Y3-D12*X2*Y 1+2*D12*X2*Y2-D12*X2*Y3+ 
. D12*X3*Y 1-D12*X3*Y2+D22*X1*X2-D22*X1*X3-D22*X2**2+ 

. D22*X2*X3))/(3*DJ) 

SKE(4, 1)=SKE( 1 ,4) 

SKE(4,2)=SKE(2,4) 

SKE(4,3)=SKE(3 ,4) 

SKE(4,4)=(4 1|t (D 1 1 * Y 1 * *2-D 1 1 * Y 1 * Y 2-D 1 1 * Y 1*Y3+D1 1*Y2**2- 

D11*Y2*Y3+D11*Y3**2-2*D12*X1*Y1+D12*X1*Y2+D12*X1*Y3+ 

. D12*X2*Y 1-2*D12*X2*Y2+D12*X2*Y3+D12*X3*Y1+D12*X3*Y2- 
. 2*D12*X3*Y3+D22*X1**2-D22*X1*X2-D22*X1*X3+D22*X2**2- 
. D22*X2*X3+D22*X3**2))/(3*DJ) 

SKE(4,5)=(4*(D 1 1*Y 1 *Y2-D1 1*Y 1*Y3-D1 1*Y2**2+D1 1*Y2*Y3- 
. D12*X1*Y2+D12*X1*Y3-D12*X2*Y 1+2*D12*X2*Y2-D12*X2*Y3+ 

. D12*X3*Y 1-D12*X3*Y2+D22*X1*X2-D22*X1*X3-D22*X2**2+ 

. D22*X2*X3))/(3*DJ) 

SKE(4,6)=-(4*(D1 1*Y1**2-D1 1*Y1*Y2-D1 1*Y1*Y3+D11*Y2*Y3 
. -2*D12*X1*Y 1+D12*X1*Y2+D12*X1*Y3+D12*X2*Y1-D12*X2*Y3 
. +D12*X3*Y1-D12*X3*Y2+D22*X1**2-D22*X1*X2-D22*X1*X3+ 

. D22*X2*X3))/(3*DJ) 

SKE(5, 1 )=SKE( 1 ,5) 

SKE(5,2)=SKE(2,5) 

SKE(5,3)=SKE(3,5) 

SKE(5,4)=SKE(4,5) 

SKE(5,5)=(4*(D1 1 * Y 1 * *2-D 1 1 * Y 1 *Y2-D 1 1 *Y 1 * Y3+D1 1* Y2**2- 
. D1 1*Y2*Y3+D1 1*Y3**2-2*D12*X1*Y1+D12*X1*Y2+D12*X1*Y3+ 

. D12*X2*Y 1-2*D12*X2*Y2+D12*X2*Y3+D12*X3*Y1+D12*X3*Y2- 
. 2*D12*X3*Y3+D22*X1**2-D22*X1*X2-D22*X1*X3+D22*X2**2- 
. D22*X2*X3+D22*X3**2))/(3*DJ) 

SKE(5,6)=-(4*(D 1 1 *Y 1 *Y 2-D1 1* Y 1 *Y3-D1 1* Y2*Y3+D1 1 *Y3**2 
. -D12*X1*Y2+D12*X1*Y3-D12*X2*Y1+D12*X2*Y3+D12*X3*Y1+ 

. D12*X3*Y2-2*D12*X3*Y3+D22*X1*X2-D22*X1*X3-D22*X2*X3+ 

. D22*X3**2))/(3*DJ) 

SKE(6, 1 )=SKE( 1 ,6) 

SKE(6,2)=SKE(2,6) 

SKE(6,3)=SKE(3,6) 

SKE(6,4)=SKE(4,6) 

SKE(6,5)=SKE(5,6) 

SKE(6,6)=(4*(D1 1*Y 1**2-D1 1*Y1*Y2-D1 1*Y 1*Y3+D1 1*Y2**2- 

. D11*Y2*Y3+D11*Y3**2-2*D12*X1*Y1+D12*X1*Y2+D12*X1*Y3+ 

. D12*X2*Y 1-2*D12*X2*Y2+D12*X2*Y3+D12*X3*Y1+D12*X3*Y2- 

. 2*D12*X3*Y3+D22*X1**2-D22*X1*X2-D22*X1*X3+D22*X2**2- 

. D22*X2*X3+D22*X3**2))/(3*DJ) 

return 

end 
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4.6.4 REDUCE program for mass matrix 


MATRIX N(1,6),SQN(6,6),NS(6,6),DMASS(6,6),SDIF(6,6); 
N:=MAT((N1,N2,N3,N4,N5,N6)); 

% — — — - 

% Calculating the integrand. 

% 

SQN:=TP(N)*N; 

% 

% Inputting shape functions 

% . 

N4:=4*S1*S2; 

N5:=4*S2*S3; 

N6:=4*S3*S1; 

Nl:=Sl-( 1/2)*(N6+N4); 

N2:=S2-(1/2)*(N4+N5); 

N3:=S3-(1/2)*(N5+N6); 

% 

% Evaluating the integration. 

% 

LET S3=1-S1-S2; 

FOR I:=l:6 DO FOR J:=l:6 DO IF Ic^J THEN « 

A:=INT(SQN(I,J),S2);B:=SUB(S2=1-S1,A)-SUB(S2=0,A); 

C:=INT(B,S1);NS(I,J):=SUB(S1=1,C)-SUB(S1=0,C)» 

ELSE NS(I,J):=NS(J,I); 

% 

% Outputting the results. 

% 

ON FORT; 

OFF ECHO; 

OFF PERIOD; 

OUT "mass.ftn"; 

WRITE ■ SUBROUTINE MASS(Xl,Yl,X2,Y2,X3,Y3,RO,CP,DMASSr 
WRITE" IMPLICIT REAL*8(A-H,0-Z)"; 

WRITE " DI MENSION DMASS(6,6) " ; 

WRITE " AREA=-(X3*Y2-X3*Yl-X2*Y3+X2*Yl+Xl*Y3-Xl*Y2)/2"- 

DMASS:=RO*CP*2*AREA*NS; 

WRITE" RETURN"; 

WRITE" END"; 

SHUT "mass.ftn"; 

% — - ----- 

% checking the correctness of results. 

%- - 

OFF FORT; OUT "CHECK"; 

R:=FOR I:=l:6 SUM «FOR J:=l:6 SUM DMASS(I,J)»; 
SDIF:=DMASS-TP(DMASS); SHUT "CHECK"; 

BYE; 


4.6.5 Fortran results of mass matrix 

SUBROUTINE MASS(X 1 ,Y 1 > X2 > Y2,X3 > Y3 ( RO,CP.DMASS) 
IMPLICIT REAL*8(A-H,0-Z) 
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DIMENSION DMASS(6,6) 

AREA=-(X3* Y2-X3* Y I -X2* Y3+X2* Y I+X 1 *Y3-X1 *Y2)/2 

DMASS( 1 , 1)=( AREA*CP*RO)/30 

DMASS(1,2)=-(AREA*CP*RO)/180 

DMASS( 1 ,3)=-(AREA*CP*RO)/ 180 

DMASS(I,4)=0 

DMASS( l,5)=-(AREA*CP*RO)/45 
DMASS(1,6)=0 

DMASS(2,1)=-(AREA*CP*RO)/180 

DMASS(2,2)=(AREA*CP*RO)/30 

DMASS(2,3)=-(AREA*CP*RO)/180 

DMASS(2,4)=0 

DMASS(2,5)=0 

DMASS(2,6)=-(AREA*CP*RO)/45 

DMASS(3, l)=-( AREA*CP*RO)/ 180 

DMASS(3,2)=-(AREA*CP*RO)/180 

DMASS(3,3)=(AREA*CP*RO)/30 

DMASS(3,4)=-(AREA*CP*RO)/45 

DMASS(3,5)=0 

DMASS(3,6)=0 

DMASS(4,1)=0 

DMASS(4,2)=0 

DMASS(4,3)=-(AREA*CP*RO)/45 

DMASS(4,4)=(8*AREA*CP*RO)/45 

DMASS(4,5)=(4*AREA*CP*RO)/45 

DMASS(4,6)=(4*AREA*CP*RO)/45 

DMASS(5, l)=-(AREA*CP*RO)/45 

DMASS(5,2)=0 

DMASS(5,3)=0 

DMASS(5,4)=(4*AREA*CP*RO)/45 

DMASS(5,5)=(8 :)! AREA*CP*RO)/45 

DMASS(5,6)=(4*AREA*CP*RO)/45 

DMASS(6,1)=0 

DMASS(6,2)=-(AREA*CP*RO)/45 

DMASS(6,3)=0 

DMASS(6,4)=(4* AREA * CP* RO)/45 

DMASS(6,5)=(4*AREA*CP*RO)/45 

DMASS(6,6)=(8*AREA*CP*RO)/4S 

RETURN 

END 


4.6.6 Checking correctness of results 

In some cases, the results from REDUCE are lengthy as the stiffness matrix shown 
above. It is very important to find a way to check their correctness. In general, a small problem 
with a known solution is used to test the correctness of a symbolic program before application 
to actual problems. In addition, some specific checking procedures are also available in each 
individual field. They require a knowledge of the specific field. Taking the mass matrix 
problem above as an example, the summation of the entries of the mass matrix should be equal 
to unit multiplied by the accessory constants. The symmetry of mass matrix is proven by 
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subtracting the mass matrix from its transpose to get zeros for each entries. The REDUCE 
program to check correctness is appended in the program shown in subsection 4.6.4. The 
following results include two parts. The first R is the summation of each entries of mass 
matrix. The second SDIF(i.j) are the difference of mass matrix and its transpose. These 
checking results confirm the correctness of REDUCE program. 

% - 

% Checking the correctness of mass matrix by summing each entries 
% of mass matrix to make an unit multiplied by accessory constants 

% - 

R := AREA* CP* RO 

% - - 

% Checking the symmetry of mass matrix by finding the difference 
% between mass matrix and its transpose. 

% 

SDIF(1,1) := 0 
SDIF(1,2) := 0 
SDIF(1,3) := 0 
SDIF(1,4) := 0 
SDIF(1,5) := 0 
SDIF(1,6) := 0 
SDIF(2,1) := 0 
SDIF(2,2) := 0 
SDIF(2,3) := 0 
SDIF(2,4) := 0 
SDIF(2,5) := 0 
SDIF(2,6) := 0 
SDIF(3,1) := 0 
SDIF(3,2) := 0 
SDIF(3,3) := 0 
SDIF(3,4) := 0 
SDIF(3,5) := 0 
SDIF(3,6) := 0 
SDIF(4,1) =0 
SDIF(4,2) := 0 
SDIF(4,3) := 0 
SDIF(4,4) := 0 
SDIF(4,5) := 0 
SD1F(4,6) := 0 
SDIF(5,1) := 0 
SDIF(5,2) := 0 
SDIF(5,3) := 0 
SDIF(5,4) := 0 
SDIF(5,5) := 0 
SDIF(5,6) := 0 
SDIF(6,1) := 0 
SDIF(6,2) := 0 
SDIF(6,3) := 0 
SDIF(6,4) := 0 
SDIF(6,5) := 0 
SDIF(6,6) := 0 
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4.7 Closed form solution of stiffness matrix of four-node element 

4.7.1 Introduction 

Although the methodology to construct the stiffness matrix of four-node isoparametric 
quadrilateral element for the plane elasticity problem is the same as that of the triangular element 
shown in last section, the techniques are quite different from each other. In the triangular 
element, the Jacobian determinant is constant and therefore the integration is straightforward. 
However the same advantage can't be gained for the isoparametric quadrilateral element. In 
general, the determinant of the Jacobian is a function of the natural coordinates. Having the of 
Jacobian determinant in the denominator of the integrand due to the coordinate transformation 
produces a tremendous difficulty in performing integration analytically, therefore the numerical 
Gauss quadrature rule is usually introduced to solve this problem. The discussions of this 
difficulty and the introduction of Gauss quadrature rule can be found in most relevant 
literatures, such as Zienkiewicz [10], Becker & Carey & Oden [11], Cook [12], Reddy [13], 
Huebner [14], Weaver & Johnson [15]. 

The inability to perform analytic integration introduces the integration error in the finite 
element results. The following paragraphs will show that this difficulty has been overcome and 
the exact closed-form solution has been obtained by appropriate application of REDUCE [8], 

4.7.2 Preliminary formulation 

The local stiffness matrix for a 2-D isoparametric quadrilateral element is formulated by 
1 1 

K - J ’ J B 1 * E *B * t *\j\d% dt) (4.45) 


Where 


• K : local stiffness matrix. 

• E : material property matrix. 

• IJI : determinant of Jacobian matrix. 

• q : natural coordinates. 

• t : thickness of element. 
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• B : strain-displacement matrix. 

In general, each entry of strain-displacement matrix B is a function of tj, IJI and can 
be expressed as : 


VS.n.UI) 


c 0 +cg +c 2 T] +c 3 | rj 

m 


(4.46) 


Where 


• bjj : are denoted by the entry in ith row and ’jth column of matrix B 

• Cj : are constants. 


For simplicity, the material property matrix E and the thickness t are assumed to be 
independent of natural coordinates. The integrand in equation (4.45) therefore will be function 
of t] and IJI, too. The entry of integrand can be expressed specifically as : 

,, d o + d £ + <*,»? +d £ 2 + dgn +d sn 1 + d£ 1 r } +d 7 & 1 2 + d£ 2 r 1 2 _ 

*„(§.»». Ub- [jj^ < 4 - 47 ) 


Where di, do, dg are constants, too. 


The Jacobian J is 
r dx 



dx 

drj 


dy ■ 




d% 


'u 

Jn 

dy 


/« 

3 11 

dr] 





(4.48) 


And the global coordinate variables x, y can be transformed to local coordinate by the same 
shape functions as those of field variables. This is an intrinsic property of an isoparametric 
element. 

xm ^ N i x i * y (4 - 49) 


Where 

• Xj, y^ : are global node coordinates. 
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• N| : shape functions. 

Specifically, the individual shape functions are 


^i-tU -£)(!-»?) (4.50) 

+ 1)0 - n) (4.5i) 

^3- + DU + V) (4.52) 

^4 ■ “ l)(l + V) (4.53) 

and the entries of Jacobian are derived from equations (4.48) to (4.53). 
j dx n, . i . 

J “~df~4 (x i~ X 2 + x 3- x J + 4(-x l +x 2 + x 3 -x 4 )-a xV +b x (4.54) 

/ ** *, 1 

J n~-^-^4< y l -y 1 + y 3 -y 4 ) + 7i-y i + y 2 + y 3 -y 4 )~a y r 1 +b y ( 4 . 55 ) 

, dx_ g i 

J u m dr) “ 4 '*i ■ K i + X 3~ x J + 'i(~ x i- x 2 +x 3 + x 4 ) - QxM + c x (4.56) 


r dX S, . i , 

J 22-~Mr-^<y l -y2 + y 3 -yJ + 4(-y l -y 2 + y 3 +y 4 )-a y z+c y ( 4 . 57 ) 
Therefore, the determinant of the Jacobian will be 

PI- 

- (b y a y -a y b y )% + (a,c y - a y c y )>) +(b t c y * b y c x ) (4.58) 

- Of + Vr) + W 

where U, V, W are independent of natural coordinates. 

Obviously, the determinant of the Jacobian is only a linear function of natural 
coordinates. This linearity allows the exact integration to be performed and the logarithm 
function is expected in the solution. As the first integration with respect to x is done, there is no 
longer an integration variables h in the denominator. Therefore the second integration with 
respect to h can also be performed analytically. However, the algebraic operations to finish 
these two integrations are too tedious to be handled by hand. Fortunately, with the help of the 
symbolic and algebraic manipulator REDUCE, these mathematical operations can be done by 
computer. In addition, the solutions can be organized in a systematic way and converted into a 
FORTRAN-code subroutine to be called by the main program. All of these procedures and 
parts of solutions are demonstrated in the next paragraphs by an example of linear elasticity. 
The explanations of commands and the time consumed in each individual command are also 
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commented for reference. The total time consumed in this execution by REDUCE is around 
two and half hours. The resultant fortran expressions for just an element of stiffness matrix 
occupy almost sixteen pages. These huge expressions are the reason why the closed form 
solution was not available before. 


4.7.3 REDUCE program and explanation 


% 

% Turning on the CPU elapse time. 

% 

1: on time; 

Time: 133 ms 

% 

% Inputting 4 shape functions. 

% p and q are natural coordinates. 
% p : xi 
%q : eta 

% - — 

2: sl:=(l-p)*(l-q)/4; 

P*Q - P - Q + 1 

51 := — 

4 

Time: 600 ms 
3: s2:=(l+p)*(l-q)/4; 

P*Q- P + Q- I 

52 := - 

4 

Time: 167 ms 


4: s3:=(l+p)*(l+q)/4; 

P*Q + P + Q + 1 

53 := 

4 

Time: 167 ms 
5: s4:=(l-p)*(l+q)/4; 

P*Q + P-Q- 1 

54 := 

4 


Time: 166 ms 


% Expressing x & y by shape function 
% and global node coordinates. 

% 
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6: x:=sl*xl+s2*x2+s3*x3+s4*x4; 

X := (P*Q*X1-P*Q*X2+P*Q*X3-P*Q*X4-P*X1+P*X2+P*X3-P*X4- 
Q*Xl-Q*X2+Q*X3+Q*X4+Xl+X2+X3+X4)/4 
Time: 400 ms 

7: y:=sl*yl+s2*y2+s3*y3+s4*y4; 

Y := (P*Q*Y1-P*Q*Y2+P*Q*Y3-P*Q*Y4-P*Y1+P*Y2+P*Y3-P*Y4- 
Q*Y 1 -Q* Y 2+Q* Y3+Q* Y 4+Y l+Y2+Y3+Y4)/4 
Time: 267 ms 




% Declaring and Inputting matrix elements to calculate the strain-displacement matrix B. 
% c : coefficient matrix 
% jac : combination of Jacobian matrix 
% sd : matrix containing the derivative of shape function. 

% b : strain-displacement matrix. 

% detj : determinant of Jacobian. 

% j I l,j 12,j21,j22 : element of Jacobian matrix. 

% 

8: matrix c(3,4),jac(4,4),sd(4,8),b(3,8)$ 

Time: 550 ms 

9: c: =mat(( 1 ,0,0,0), (0,0,0, 1 ) ,(0, 1 , 1 ,0))$ 

Time: 183 ms 

1 0: jac: =mat((j22,-j 1 2,0,0), (-j2 1 ,j 1 1 ,0,0) , 

(0,0,j22,-j 12),(0,0,-j2 1 j 1 1))$ 

Time: 284 ms 

1 1: sd:=mat((df(sl,p),0,df(s2,p),0,df(s3,p),0,df(s4,p),0), 
(df(sl,q),0,df(s2,q),0,df(s3,q),0,df(s4,q),0), 
(0,df(sl,p),0,df(s2,p),0,df(s3,p),0,df(s4,p)), 
(0,df(sl,q),0,df(s2,q),0,df(s3,q),0,df(s4,q)))$ 

Time: 833 ms 

12: b:=c*jac*sd/detj$ 

Time: 417 ms 

% 

% Inputting material matrix D and calculating integrand. 

% D : material matrix (assumed symmetric). 

% Ga : integrand. 

% th : thickness of element. 

%- 

13: matrix d(3,3),ga(8,8); 

Time: 316 ms 

14: d:=mat((ell,el2,el3),(el2,e22,e23),(el3,e23,e33))$ 

Time: 200 ms 

15: ga:=tp(b)*d*b*th*detj$ 

Time: 8067 ms 
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% - 

% Evaluating each element of Jacobian. 

% 

16: on factor; 

Time: 150 ms 

17: on div; 

Time: 33 ms 

18: on rat; 

Time: 50 ms 

19: jl l:=df(x,p); 

1 

Jl 1 := — *((X1 - X2 + X3 - X4)*Q - XI + X2 + X3 - X4) 

4 

Time: 750 ms 

20: jl2:=df(y,p); 

1 

J12 := — *((Y 1 - Y2 + Y3 - Y4)*Q - Y1 + Y2 + Y3 - Y4) 

4 

Time: 550 ms 

21: j21:=df(x,q); 

1 

J21 := - — *((X1 + X2 - X3 - X4) - (XI - X2 + X3 - X4)*P) 
4 

Time: 667 ms 

22: j22:=df(y,q); 

1 

J22 := - — *((Y 1 + Y2 - Y3 - Y4) - (Y 1 - Y2 + Y3 - Y4)*P) 

4 

Time: 650 ms 


% 

%Making substitution for Jacobian matrix. 

% ax=(xl-x2+x3-x4)/4, bx=(-xl+x2+x3-x4)/4 
% ay=(yl-y2+y3-y4)/4, by=(-yl+y2+y3-y4)/4 
% cx=(-xl-x2+x3+x4)/4, cy=(-yl-y2+y3+y4)/4 

% 

23: jll:=ax*q+bx; 

Jll := AX*Q + BX 
Time: 333 ms 

24: jl2:=ay*q+by; 

J12 := AY*Q + BY 
Time: 1 17 ms 

25: j21:=ax*p+cx; 

J21 := AX*P + CX 
Time: 1 17 ms 
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26: j22:=ay*p+cy; 

J22 := AY*P + CY 
Time: 116 ms 

% 

%Calculating the determinant of Jacobian. 

%— 

27: matrix j(2,2),ske(8, 8); 

Time: 200 ms 

28:j:=mat((jll,jl2),(j21j22)); 

J(l,l) := AX*Q + BX 
J( 1,2) := AY*Q + BY 
J(2,l) := AX*P + CX 
J(2,2) := AY*P + CY 
Time: 367 ms 

29: detj:=det(j); 

DETJ := - ((AX*P + CX)*(AY*Q + BY) - (AX*Q + BX)*(AY*P + CY)) 
Time: 333 ms 


% 

%Making a further substitution and giving the lineality of determinant. 

% 

30: let -ax*cy+ay*cx=al,ax*by-ay*bx=a2,-bx*cy+by*cx=a3; 

Time: 450 ms 

31: detj:=detj; 

DETJ := - (A1*Q + A2*P + A3) 

Time: 450 ms 

32: on exp; 

Time: 50 ms 

%--- 

%Performing the double integration. 

%Due to the complication of the expression in the integrand, it is necessary to make the "pre- 
%treatment" to each element of integrand before integration. This is a vital step to avoid the 
^limitation of memory space and finish the job. 

33: for i:=l:8do forj:=l:8do 
if j>=i then 

«cp:=coeff(ga(i,j),p); 
p0:=(first cp); 
pl:=( second cp); 
p2:=( third cp); 
low:=den(ga(i,j)); 
ga(i,j):=(dOfdl*p4-d2*p**2)/low; 
cl:=int(ga(ij),p); 
c2:=sub(p= 1 ,cl )-sub(p=:- 1 ,c 1) ; 
c3 : =sub{d0=p0* low.d 1 =p 1 *low,d2=p2*low,c2) ; 
c4:=sub(log(al *q-a2+a3)=mlog,log(al *q-fa2+a3)=plog,c3); 

on exp; 
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cq:=coeff(c4,q); 
qO:=(first cq); 
ql:=( second cq); 
q2:=( third cq); 

kO 1 :=lcof(num(qO),mlog)/den(qO); 
qO:=reduct(num(qO),mlog)/den(qO); 
k02:=lcof(num(q0),plog)/den(q0); 
k03 : =reduct(num(qO),p!og)/den(qO) ; 
k04:=lcof(num(ql),mlog)/den(ql); 
q 1 :=reduct(num(q 1 ) ,mlog)/den(q 1 ) ; 
k05:=lcof(num(q 1 ),plog)/den(q 1 ) ; 
k06:=reduct(num(q 1 ),plog)/den(q 1 ) ; 
k07:=lcof(num(q2),mlog)/den(q2); 
q2:=reduct(num(q2),mlog)/den(q2); 
k08:=lcof(num(q2),plog)/den(q2); 
k09:=reduct(num(q2),plog)/den(q2); 
c4: = dO 1 * 1 og(a 1 * q -a2+a3 ) +d02 * I og(a 1 * q+a2+a3 )+d03 
+d04*q*log(a 1 *q-a2+a3)+d05*q*log(a 1 *q+a2+a3)+d06*q 
+d07*q**2*log(al*q-a2+a3)+d08*q**2*log(al*q+a2+a3) 
+d09*q**2; 
c5:=int(c4,q); 

c6: =sub(q= 1 ,c5)-sub(q=- 1 ,c5) ; 

let log(al**2-a 1 *a2+a 1 *a3)=h 1 ,log(-al * *2-al *a2+al *a3)=h2, 
1 og(a 1 +a2+a3 )= h3 , 1 og( -a 1 +a2+a3 ) =h4, 1 og( a 1 -a2+a3)=h5, 
log(-a 1 -a2+a3)=h6; 
factor hl,h2,h3,h4,h5,h6; 

ske( i ,j ) : =sub(dO 1 =k0 1 ,d02=k02,d03=k03 ,d04=k04,d05=k05 
,d06=k06 ) d07=k07,d08=k08,d09=k09,c6)» 
else ske(i,j):=ske(j,i); 

Time: 234766 ms 


% 

%Showing the results for the element in 1st row and 1st column of the local stiffness matrix. 
%Hl=log(al**2-al*a2+al*a3), H2=log(-al**2-al*a2+al*a3) 

%H3=log(al+a2+a3) , H4=log(-al+a2+a3) 

%H5=log(al-a2+a3) , H6=log(-al-a2+a3) 

% - 

34: ske:=ske; 


1-1 1-1-12 
SKE(1,1) := H1*TH*(— *A1 *A2*E13-— *A1 *A2 *A3 *E13 + 

8 8 

1 - 1-1 2 1-1 -1 

— *A1 *A2 *A3*BX *E33 *A1 *A2 *A3*BX*BY*E13 

16 8 

1-2 2 1 

+ -— *A2 *A3 *E13 + — -*E13) 

16 16 

12-2 12-32 

+ H6*TH*( — *A 1 *A2 *E13 + — *A1 *A2 *AX *E33 - 
16 48 
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1 l -1 1 -1 

+ — - -*E13) +TH*( - *A1*A2 *E13 + — *A1 *A2*E13 + 

16 4 4 

1-2 1-2 2 

-— *A2 *A3*BX*BY*E13 + — *A2 *A3*BY *E11) 

3 6 

Time: 79317 ms 


4,7.4 Fortran subroutine from REDUCE 

% 

%Converting the results into FORTRAN code 

% 

subroutine(xl ,y 1 ,x2,y2,x3,y3,x4,y4,el t e2,e3,e4,e5,e6,th,ske) 
implicit real*8(a-h,o-z) 
dimension ske(8,8) 
ax=(xl-\2+x3-x4)/4, 
bx=(-x 1 +x2+x3-x 4)/4. 
a y=(yl-y2+y3-y4)/4. 
by=(-yl+y2+y3-y4)/4. 
cx=(-xl-x2+x3+x4)/4. 
cy=(-yl-y2+y3+y4)/4. 
al=-ax*cy+ay*cx 
a2=ax*by-ay*bx 
a3=-bx*cy+by*cx 
h 1 =log(a 1 * * 2-a 1 *a2+a 1 * a3) 
h2=log(-al**2-al*a2+al*a3) 
h3=log(al+a2+a3) 
h4=log( -a 1 +a2+a3) 
h5=log(al-a2+a3) 
h6=log(-al-a2+a3) 

ANS 14= 1/4* A 1 **(-3)*A3**2* AY *CX*E13- 1/8* A 1**(-3)*A3**2 
. *AY*CY*Ell-l/16*Al**(-3)*A3**2*CX**2*E33+l/8*Al**(-3 
)*A3**2*CX*CY*E13-1/16*A1**(-3)*A3**2*CY**2*E1 1-1/16 
. *A2**(-2)*A3**2*E13+1/16*E13 


ANS1=H1*TH*ANS2 

ANS28=-1/4*A1**(-3)*A3**2*AY*CX*E13+1/8*A1**(-3)*A3** 
. 2*AY*CY*E1 1 + 1/16*A1**(-3)*A3**2*CX**2*E33-1/8*A1**( 
-3)*A3**2*CX*CY*E13+1/16*A1**(-3)*A3**2*CY**2*E1 1+1/ 
16*A2**(-2)*A3**2*E13-1/16*E13 


ANS 15=H2*TH* ANS 16 

ANS48=1/8*A2**(-3)*A3**2*BX*BY*E13-1/16*A2**(-3)*A3** 

2*BY**2*E11 


ANS29=H3*TH*ANS30 

ANS68=-1/16*A2**(-3)*A3**2*BX**2*E33+1/8*A2**(-3)*A3 
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** 2 *bX*BY*E13- 1/16*A2**(-3)*A3**2*BY**2*E1 1 


ANS49=H4*TH*ANS50 

ANS75=1/16*A2**(-3)*A3**2*AX**2*E33-1/8*A2**(-3)*A3** 

. 2*AX*AY*E13+1/8*A2**(-3)*A3**2*AX*BX*E33+1/16*A2**( 
-3)*A3**2*AY**2*E1 1- l/4*A2**(-3)*A3**2*AY*BX*E13+l/8 
. *A2**(-3)*A3**2*AY*BY*E1 1 + 1/16*A2**(-3)*A3**2*BX**2* 
. E33-1/8*A2**(-3)*A3**2*BX*BY*E13+1/16*A2**(-3)*A3**2 
. *BY**2*E11-1/16*E13 


ANS69=H5*TH*ANS70 

ANS82=-3/16*A2**(-2)*A3*BY*CY*Ell+l/16*A2**(-3)*A3**2 
*AX**2*E33- l/8*A2**(-3)*A3**2*AX*AY*E13+l/8*A2**(-3) 
*A3**2*AX*BX*E33+1/16*A2**(-3)*A3**2*AY**2*E1 1-1/4* 
A2**(-3)*A3**2*AY*BX*E13+1/8*A2**(-3)*A3**2*AY*BY* 

. El 1+1/16* A2**(-3)*A3**2*BX**2*E33-1/8*A2**(-3)*A3**2 
. *BX*BY*E13+1/16*A2**(-3)*A3**2*BY**2*E1 1+1/16*E13 


ANS83=TH*ANS84 

ske( 1 , 1)=ANS 1+ANS 1 5+ANS29+ANS49+ANS69+ANS76+ANS83 


return 

end 


4.8 Significance and conclusion 

Some conclusions are drawn and the significance of automatic problem formulation is 
discussed as follows : 

1. Improving on-line efficiency — the closed-form solution of the local stiffness matrix 
allows us to get a numerical value by simply substituting the global nodal coordinates 
into a FORTRAN-code subroutine. This procedure is done in just one step. Of course, 
this is faster than the Gauss quadrature rule which usually needs more than one 
integration point to get a reasonable solution 4 . The symbolic template in the nonlinear 
numerical analysis also plays the same role. In the case of a large number of elements (or 
large dimension size in matrix), the significance in improving on-line efficiency will be 
greater. 

2. Increasing the accuracy of solution — the closed-form solution is an exact solution. 
There is no integration error introduced into the evaluation of the stiffness matrix. The 

4 Reduce integration is an exception and sometimes results in Hourglass drawback. This 
special case is ruled out here. 
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accuracy of the Fourier series approximation can be increased as high as the user 
requires. Therefore, the results from using symbolic and algebraic manipulation will be 
more precise than those of pure numerical analysis. 

3. Free from hand-calculation and typing error — As the results in the tensor formulation, 
derivation of equations of motion and stiffness matrix constructions show, the algebraic 
expressions are too lengthy to be formulated by hand. Even supposing that they could be 
done by hand, it would be so tedious that nobody could guarantee that no mistakes 
would be made when trying to key them into the computer. With the use of symbolic 
and algebraic manipulation, both difficulties are completely solved. As long as the user 
inputs the correct commands, there will not be any question about the correctness of the 
results. 

4. Simplifying FORTRAN programming — the numerical values of the local stiffness 
matrix can be obtained by simply using the "CALL" command once. This isn't true 
when Gauss quadrature integration is employed in the finite element method to evaluate 
integration. It is necessary to have a "DO" loop, "CALL" command and multiplications 
of various weight coefficients for different integration points. These will complicate the 
program and therefore will produce more error sources. 

5. Further analysis of symbolic results becomes available — Sometimes the pre-analysis of 
the expressions produced from symbolic and algebraic manipulation will lead to a 
dramatic improvement in the incoming numerical analysis. The closed form solution 
makes this analysis feasible. For example, suppose that the diagonal terms of global 
stiffness matrix need to be more dominant to improve the ill-condition, this can be 
achieved by appropriately relocating nodes so that the off-diagonal terms of local 
stiffness matrix will be smaller or even vanished. This is the unique advantage that the 
numerical method does not possess. 
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CHAPTER V 


APPLICATION OF SYMBOLIC AND ALGEBRAIC MANIPULATION TO A 
MATERIALLY NONLINEAR PROBLEM — RIGID-PLASTIC RING 

COMPRESSION 


5.1 Introduction 

The application of the finite element method to the rigid plastic problem was originally 
devised by Lee and Kobayashi ([1], [2]), and became popular in the last decade. This method 
allowed the deformation behavior of metal to be revealed on the computer screen stage by stage 
during the process of metal forming. As a consequence, the design techniques of die, and the 
manufacturing process were improved. This contribution to the industrial manufacturing field 
is recognized to be very significant. 

Starting from the principle of virtual work and associating with the normality condition 
of plasticity, the theoretical analysis of this method leads to an inequality objective function 
with an equality constraint. This equality constrained problem is then changed into an 
unconstrained problem by introducing Lagrange multipliers. As the stationary requirement is 
reached, the total unconstrained problem can be solved incrementally by the finite element 
method and the upper bound solution will satisfy the equilibrium equations, constitutive 
equations, compatibility equations, incompressibility constraint, and boundary conditions. 

Despite the success of the finite element application to the problem, the formulation of it 
is very tedious. Especially when the friction boundary condition is considered, the performance 
of integration along the interface surface and the evaluation of the first and second derivatives 
with respect to velocity fields is always difficult to obtain by hand. Therefore, unlike the 
traditional hand derivation, this chapter will utilize the symbolic and algebraic manipulator 
REDUCE to do the job of formulation. The quantities involved in the formulae then can be 
made into subroutines symbolically at the the element level to facilitate the global assemblage. 
As a result, otherwise intractable tasks become possible and free of errors. 
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5.2 Preliminary formulation 

Consider a body of volume V with the essential boundary condition , velocity £f 
prescribed on surface S u and two separated natural boundary conditions , traction f and 
frictional stress f, prescribed on Sp and S c , respectively [Figure 5.1]. The actual stress and 
velocity fields will satisfy the following relationships : 

1. Equilibrium conditions : 


°, t , “0 (neglecting body force ) 



Figure 5.1: Configuration of domain and boundary conditions 

2. Compatibility and incompressibility conditions: 

- x(u. + u ) 
e -0 

tt 

3. Stress-strain rate relationship 



(5.1) 


(5.2) 

(5.3) 


(5.4) 


Where 

• a ,J : deviatoric stress. 

*° and £ : effective stress and strain-rate, respectively. 

4. Boundary conditions : 
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Where the relative velocity between die and deforming body is defined as follows : 

V r - <r - - v r t~ 

The t here is the unit base vector along die and working piece interface. 

With an admissible velocity field & , the virtual work principle gives 
f v °ij E ,j dv “ f s F • udS + f f • (V r + u*)dS + f (a n )U dS 

F * e ** S u '* 1 i 

The frictional stress f* is defined as follows and is plotted in Figure 5.2. 


f (V) 

r 


i 

, 

mg 


-mg 

1 



V 


Figure 5.2 : Plot of frictional stress vs. relative velocity 

FT 


y 

fPr) - - m * g * T-r-rr 


(5.5) 

(5.6) 

(5.7) 


(5.8) 


(5.9) 


(5.10) 


where the friction factor m is in the range of 0 <; m s 1 and g is the yield shear stress. By 
considering the normality condition of yield surface 


(ct* - a )i\ a 0 

V 17 ' i; 


(5.11) 
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and inequality equation 5 


p'-mv p ; i0 

The equation (5.9) becomes 


(5.12) 


l d ', dv - S, e ■ - f A O • 

* 1 s}°‘J n < W i dS + / f ( V r ) - a* ds ( 5 . 13 ) 

* 

All conditions will be met when the left hand side of equatton (5.13) reaches the 
mtntmum values wtth respect to t and sattsfies the incompressibility condition of equation 

( d h, r r~' S "° an,b ' gU ' ly omittin « ,he “"risk, for the sake of simplicity the 
mtsstble field u wtll be denoted as u for the follow, „g dtscusstons without any special 

"mo ,he y o n b‘T m f 8 ‘ he U8ra " 8e mUU ' P " er '• ‘ he K ' uat '™ *» b. included 

therefore ° n ' problera f ° r eleme " 1 formulation * 

" aF l -C S * V - I/’ *« ~ //< ?,)•?,« + JT Xe, dV | . o (5.14) 

Iff re V k 7 PreSentS lhe flmC,ional inside the square bracket Since the frictional stress is not 
t erentta e at thepomt V, -0 [see Figure 5.21, dir" does not exist and the convergent 

d";°- ; Va " able f ° r equa,i0n (514 >- •» -«• •» overcome thts numerical 
Ity. fncttonal stress is approximated in terms of arctangent functton [see Figute 5.3) 


f (V r ) =* - m *g *[(^-)ton - X 


r 

a\i/ 


-)]/' 


(5.15) 


where the arbttrary constant a is several orders less rhan dm die velocity. I. functton is to 
exaggerate the argument of arctangent to reduce the error between equation ,5., 5) and equatton 
• 0).Theequauon (5. 15) is absorbed into functional by the follow, ng inequality* 


^ See appendix A for proof. 
See appendix B for proof. 
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Figure 5.3 : Arctangent approximation of frictional stress 

jj ^ ' dV ' ~f 0 *f-dV r *f{V,)-(V' r -V T ) 

The final form of equation (5. 13) is 

^ d f* f* J p* I 

~K ■ ifl f dtdV - J,f -fids- J s lf YdV,)dS * fxe.dV ] . 


z 



(5. 16) 


0 (5.17) 
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5.3 Matrix method for the ring compression problem 


If a specific problem of ring compression is considered [see Figure 5.4], the surface 
traction is due to friction force only. Therefore the term for S f integration in equation (5. 17) 
can be dropped. By the substitution of the following equations to equation (5. 17), 

{“}-£*£/ (5.18) 


(5.19) 

the stationary requirement ot the tunctional results in the following nonlinear equation. 

/ (-^F 2 -)* *ZdV + JV *K *U *^z^dV + A *Q 

^dV,)dS ( 5 . 20 ) 

together with an incompressibility constraint equation 

V T *Q - 0 (5.21) 


'9 

Yu 


-VI 


± 

dr 


0 


0 


r 0 


i 

r 

d_ 

.dz 


e e 


dz 


d_ 

dr 




where 

• K - B T *D *B 

,Q - fB T *C dV 

v 

T • 

• C : proper matrix such that implies the C e incompressibility condition. 

• D : flow matrix. 

•U: vector of velocity. It is represented as & previously. 

The derivative part in the second term of equation (5.20) is equal to 
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(5.22) 


^ e l 1 do a dt Y de 

dU “ T dU ft dU ~ ~ -p~dU 

The first term of (5.22) is dropped for simplicity due to the assumption that a is constant for 
the infinitesimal strain change [see References 2], Physically, a is the current yield strength 
(denoted as Y) of the material. 


The equations (5.20) and (5.21) are then perturbed by introducing a small velocity 
Af/ into the velocity vector. This gives the following equations : 


I v + ^-(-y-)A U ]*K *(U + A U)dV + A * Q 

+ f(U + AU) T * K *(U + AU )[-£-(%.) + JL^^SL^u ]dv 

dU s dU c 


r d® 

[ dU + 


dJP 

dU 1 


AU] 


(5.23) 


and 


(U + AU j *Q = 0 


where 


(5.24) 


® 


-/</' 


F dV r )dS 


is contributed by the friction traction. After neglecting the second and higher order terms, the 
equations (5.23) and (5.24) can be combined into the following forms : 


( \P Q 


d 2 ® 1 


AU 

q t 0 

+ 

dU 

4 

A > 

\ L J 


0 0 

/ 

Y 



_1 

Y 


d®_ 
dU . 
0 


J 9.1 


(5.25) 


where 

P ,;- 2 l{j K -<4r)‘/n} 

h ™-S.{t k ' u -j e } 
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(5.28) 

(5.29) 


2 Hi Ti 

£ “ E r + £ z + 


e 

E 


Eg E r e z £ z £ g £ r £ g + 77 r Z 


8*1 


B T * A 


e,„- fB'cdv 


A - r i j_\ 

4 *i l Se ‘ 4 ’ 1 9f ' 


2£ r " - £ fl 

2l! ‘< - £ 9 - £ r 

2£ ' 8 " £ r - ** 
3 • 

TT,, 


* D • B 

M,„ - < 7 ) ! [ j -£ -K *V ]L r 


A - ;C r if 

- £ *f/ r ♦* 


k -U T * K *U 

_ f J'l 2 . V r 

0 - - f s [f Q m*g* r ]dS 

m, a g are constants for a specific material. 


5.4 Finite element analysis 


(5.30) 

(5.31) 


(5.32) 

(5.33) 

(5.34) 

(5.35) 

(5.36) 

(5.37) 


The domain discretization for ring cross section is shown in Figure 5.5. The actual 
domain used by finite element analysis is only the upper half of that given in Figure 5.5 due to 
the symmetric geometry. 


The finite element model contains 96 four-node quadrilateral elements with 1 17 nodes 
in total. For an isoparametric element, the shape functions are 


?1 -1(1-5)(1-/) 

tf 3 “ 4<1 +*)(l + 0 


, q 4 - j(l-s)(l+0 


and 
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I COORDINATE (1N> 


Where 


u(s,t) - 2 q , * «, 

i-i 

4 

r(s,t )-^q, *r t 

i -1 


, v(s,t) - 2<7, * v , 

l -1 

4 

. z(s,t) = 2^. * z , 

i -1 


• s, t are natural coordinates. 

• r, z are physical coordinates. 

• u, v are the velocity components in the r, z direction, respectively. 

The subscripts in r, z, u, v represent the nodal indices. 

The strain rate in axisymmetric case can be expressed in a matrix form as follows 



( 1.0 1.0 2.0 3.8 


R COMOINMt <IN> 


Figure 5.5 : Discretization of ring cross section. Due to the symmetry 
of geometry, only the upper half of this cross section is 
used for finite element analysis. 
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Where 


B 


r 

kj 


rj_ 

dr 

0 

J_ 

r 

_d_ 

dz 


0 

d 

dz 

0 

£_ 

dr 


{:} 


B *U 


(5.38) 
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0 

i 
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0 

0 
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(5.39) 
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0 
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ds 
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0 

9., 

0 
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9 3l , 
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9„ 
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0 

9 3 , 
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(5.40) 


The Jacobian is 


The flow matrix is 



dr 


dz 



J - 

ds 


ds 




dr 


dz 




dt 


dt 




2 

0 





3 


0 

0 

D - 

0 

2 

3 


0 

0 


0 

0 


2 

r 

0 


0 

0 


0 

1 

1 

C r - 

(1 

1 


i 

0) 


(5.41) 


(5.42) 


(5.43) 
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Since bilinear velocity distributions are assumed for four-node quadrilateral element, 
the relative velocities along the interface are interpolated as follows : 





(5.44) 


If die velocity is specified as u ldl =l, then the frictional stress will be 


F --{m *g * (£)tan _, (-^-)}f 


(5.45) 


and 


<*>=-/ [jm *g *(r)tan r ]dS 


(5.46) 


5.5 Application of symbolic manipulation 

The difficulty of formulating the equations shown in the last paragraph can be eased by 
the employment of symbolic and algebraic manipulation. Based on equation (5.25), the original 
goal is to make a template form in element level for global assemblage. However, due to the 
limitation of memory capacity in hardware systems, the goal is modified to make a template 
form for individual entries of equation (5.25) only. There are three kinds of forms produced by 
REDUCE: 


1. Integrable form — This is the form which results from the fact that integration can be 
performed analytically by REDUCE. The evaluation of equations (5.30) and (5.37) 
belongs to this class. There are no natural coordinates in the resultant expressions. 

2. Non-integrable form — The equations (5.26), (5.27) and (5.33), for example, are not 
integrable due to the existence of £ , its square and even cubic in the denominator of the 
integrands. Moreover, since the resultant expression of £ occupies more than sixteen 
pages, the complete integrands of equations (5.26) and (5.27) are not available. 
Therefore, the individual form for quantities, such as K, M, N, k and E, are obtained 
symbolically, and the summation as well as their integration are carried out numerically 

3. Miscellaneous forms — Other equations which have nothing to do with integration, such 
as (5.28), (5.29), (5.31), (5.32), (5.35) and (5.36), are obtained symbolically by matrix 
operations. 
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The evaluation of ^ 7 - in equation (5.25) requires a number of delicate pre-treatments. It 
is necessary to discuss this subject independently here so that the fundamental nature of 
symbolic and algebraic manipulation will be revealed. Intuitively, there are three steps to 
evaluate They are : 

1. Evaluation of the integral with respect to relative velocity V f \n (5.37). 

2. Evaluation of the integral with respect to interface domain S in (5.37). 

3. The results then are differentiated with respect to velocity fields. 


Theoretically, there is nothing wrong with the methodology given above. In fact, 
REDUCE can only do the first evaluation. The other two are not feasible. Therefore, some pre- 
treatments are necessary to make REDUCE work to evaluate the This will force us to 
deviate from the formal methodology given above as follows. 

d<P d c r v ' 1 . v 

■sr-’srS.tf. ■*** v ')" 

~-m *g *(hS 'irif tan ~\ V -r)dV , 

' 0 I 

- - m * g *(x)/ s [tan-‘(^-)-^-]<i5 (5.47) 

The term — can be evaluated from equation (5.44). Equation (5.47) becomes 

d<p 2 r 2 * r ' 1 . v, <• -r 

^ - - m *g * (Z-)J o J r tan ->(-f > 7777 rdrdd 

- \r* tan^T-) - r r 2 tan'^)** (5.48) 

1 


Rewrite equation (5.44) into 


V r 



Yr +Z 


(5.49) 


and substitute it in equation (5.48). 


J<P 

du x 


- 4 * m* g 




2 . -i,r z . 
tan (— r + t) -r r. 


tan ‘‘(rr + i-)]dr 


(5.50) 


Since REDUCE is still unable to handle equation (5.50), further treatments are required. Let 
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Ti/ y z 

W - -s-r + T 


(5.51) 


Y r + z 


then 



(5.52) 


and equation (5.50) becomes 


d0 


-4* m • g 


f 

uj 


! r 0V -Z') 

[ ; — tan 


VV 


^P-^tan* 1 VP ]d\V 


(5.53) 


Note that the upper and lower limits of integration are changed. Starting from this 
point, REDUCE can proceed by itself. The tasks it performs in this specific problem include 

d<P 

the integration and back substitution of the variables. The other quantity, such as ^ etc., can 
be calculated in the similar manner. The second derivatives are simply computed by 
differentiation of the results of the first derivatives. The REDUCE program and a part of 
fortran solution are presented as follows. 

% ----- - - 

%REDUCE program to calculate friction part and its derivatives 

% - - - — 

OFF EXP; 

A1:=INT(ATAN(W)*(W-ZP)**2/YP-ATAN(W)*(W-ZP)*R2,W)*4*TM*TK/((R1- 

R2)*YP**2); 

A2:=IOT(ATAN(W)*R1*(W-ZP)-ATAN(W)*(W-ZP)**2/YP,W)*4*TM*TK/((R1- 

R2)*YP**2); 

LET YP=(ui-U2)/(A*(Rl-R2)),ZP=(Rl :<c U2-R2*Ul)/(A*(Rl-R2)); 

LET LOG((A**2+U2**2)/A**2)-LOG((A**2+Ul**2)/A**2) 
=LOG((A**2+U2**2)/(A**2+Ul**2)); 

OFF EXP; 

ON FORT; 

OFF ECHO; 

CARDNO!*:=10; 

OUT "look.ftn"; 

WRITE 

■ SUBROUTINE FRIC(U1,U2,A,TK,TM,R1,R2,B1,B2,B11,B12,B22)"; 

WRITE " IMPLICIT REAL*8(A-H,0-Z)"; 

B 1:=SUB(W=U2/A,A 1 )-SUB( W=U 1/A ,A 1) ; 

B2:=SUB(W=U2/A,A2)-SUB(W=U1/A,A2); 

B1 1:=DF(B l.U 1); 

B22:=DF(B2,U2); 

B12:=DF(B1,U2); 

WRITE" RETURN"; 

WRITE" END"; 

OFT FORT ; 

SHUT "look.ftn"; 

BYE; 
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nnn 


The fortran subroutine made by REDUCE for frictional part. 


SUBROUTINE FRIC(Ul,U2,A/rK/rMJU^31323U.B12,B22) 
IMPLICIT REAL*8(A-H.O-Z) 

ANS4^LOG((A**2+U2**2)/(A**2+Ul**2))*A**3*R2**2-3 *LOG 
. ((A**2+U2**2)/(A**2+U1**2))*A*R1**2*U2**2+3.*LCK3((A 
. **2+U2**2)/(A**2+Ul**2))*A*Rl*R2*Ul*U2+3.*LOG((A**2+ 

. U2**2)/(A**2+UI**2))*A*Rl*R2*U2**2'3.*LOG((A**2+U2** 

, 2)/(A**2+Ul**2))*A*R2**2*Ul*U2+A*Rl**2*Ul**2-6.*A*Rl 
. **2*Ul*U2+5.*A*Rl**2*U2**2+A*Rl*R2*Ui**2+6.*A*Rl*R2* 

. Ul*U2-7.*A*Rl*R2*U2**2-2.*A*R2**2*Ul**2+2 *A*R2**2* 

, U2**2 

ANS3=-6.*ATAN(U2/A)*A**2*R1**2*U2+3.*ATAN(U2/A)*A**2* 
.R1*R2*U1+9.*ATAN(U2/A)*A**2*R1*R2*U2-3,*ATAN(U2/A)*A 
. **2*R2**2*U1-3.*ATAN(U2/A)*A**2*R2**2*U2+2.*ATAN(U2/ 
A)*R1**2*U2**3-3.*ATAN(U2/A)*R1*R2*U1*U2**2-ATAN(U2/ 

. A)*R1*R2*U2**3+3.*ATAN(U2/A)*R2**2*U1*U2**2-ATAN(U2/ 

. A)*R2**2*U2**3+LOG((A**2+U2**2)/(A**2+Ul**2))*A**3* 
Rl**2-2.*LOG((A**2+U2**2)/(A**2+UI**2))*A**3*Rl*R2+ 

. ANS4 

ANS2=6.*ATAN(U1/A)*A**2*R1**2*U2-3.*ATAN(U1/A)*A**2* 

. R1*R2*U1-9.*ATAN(UI/A)*A**2*R1*R2*U2+3.*ATAN(U1/A)*A 
. **2*R2**2*U1+3.*ATAN(U1/A)*A**2*R2**2*U2-2.*ATAN(U1/ 
A)*Rl**2*Ul**3+6.*ATAN(Ui/A)*Rl**2*Ul**2*U2-6.*ATAN( 

. U1/A)*R1**2*U1*U2**2+ATAN(U1/A)*R1*R2*U1**3-3.*ATAN( 

. U1/A)*R1*R2*U1**2*U2+6.*ATAN(U1/A)*R1*R2*U1*U2**2+ 

. ATAN(U1/A)*R2**2*U1**3-3.*ATAN(U1/A)*R2**2*U1**2*U2+ 

. ANS3 

ANS1=2.*ANS2*TK*TM 

B1=ANS1/(3.*(U1-U2)**3) 


ANS 1=2.*ANS2*TK*TM 

B12=ANS1/(U1-U2)**4 

RETURN 

END 
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5.6 Numerical evaluation and result treatment 

The numerical scheme employed is the Newton-Raphson iteration with a displacement 
increment of 0.01 in each stage. The initial guesses are slightly modified from the solution of 
the elastic ring compression problem. Since the existence of zero velocities in r-direction at 
interface nodes will overflow the subroutine FRIC, the problem is modified by assigning a 
different small number to the r-component of every relevant node. The question that arises is 
how small they should be. According to the experiments, only the numbers which are smaller 
than 10 12 will achieve convergence with this scheme. The convergence criterion used here is 

0.00005 . 

The boundary conditions are specified at two parts of the boundary : 

1. Symmetric boundary condition — The velocities in z direction [see Figure 5.5] are 
specified to be zeros along the z=0 boundary. 

2. Die velocity boundary condition — The velocities at the interface surface between the die 
and the working piece are specified as unit per second in the negative z direction. 

Numerical integration of non-integrabie terms is performed by the 4-point Guass 
quadrature rule. The assembly of a global matrix is also done numerically. The equation solver 
is the Gauss elimination method, from the 1MSL subroutine library. 

The deformed configurations for friction factor for m=0.5 and m=0.0 are shown in 
Figures 5.6 and 5.7, respectively. As the figures show, the deformed shapes are completely 
different for low and high frictional factors. The velocity distributions in the deformed states 
are also plotted in Figures 5.8 and 5.9, respectively. The neutral lines in both cases are visible 
from pictures. Figures 5.10 and 5.11 also show the effective stress^ distributions for two 
cases. As the shapes of elements are distorted, the error increases and the convergence of the 
scheme becomes harder. In order to continue the execution, the technique of adaptive mesh 
needs to be introduced. 


^ The effective stress is defined as 

l 

o t — (cr* +c 7 y + o z — o x Oy — OyO z — o x o z + 3t^ + 3t„ + 3(z xl ) 
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CHAPTER VI 


APPLICATION OF SYMBOLIC AND ALGEBRAIC MANIPULATION TO THE 

PLATE PROBLEM 


6.1 Introduction 

Although many shell theories exist in literature, there are only two distinct concepts 
from which these theories are derived. One takes the three-dimensional body as a starting point 
and tries by various means to reduce the problem to the form that can be expressed in a two 
dimensional manifold. This class of theories are called derived theories. The works by W T 
Ko,ter. and E. Reissner are in this class. The theory adopted by this study belongs to this 
class. The other class of theories consider a shell as a bounded region of some deformable two 
dimensional manifold, and are supplemented by one or more fields of vectors over this 
manifold. This class of theories are called direct theories. A. E. Green, P. M. Naghdi, and W. 
L. Wain wright worked on this class of shell theories. Since the real shell is a three-dimensional 
body, the direct approach has to rely upon some a priori statements. 

Despite a large amount of publications using the finite element method to solve plate 
and shell problems, none deals with the problems by employing the tool of symbolic and 
algebraic manipulation. This is not only due to the late availability of software, but also due to 
the capacity limitations existing in the symbolic and algebraic softwares . This chapter outlines 
the simple and universal methodology to solve the plate and shell problems, then presents 
examples which apply symbolic and algebraic manipulation to them, and finally switches to a 
numerical method at the point where the symbolic manipulation is stuck by its limitations. As a 

consequence of this work, the analytic study of plate and shell problems by computer are 
pushed a step further. 


S “geScS^Tu™ de,ai ' d,SCUSSI ° n ° n ,he CapaCi ‘ y ol > mbo " c “ d 
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6.2 Preliminary formulation 


Stating from the equations of equilibrium in three-dimensional state, 
Z^a"(« a ,z) + P'(u a ,z ) -0 

Where 


•Dj : operator of covariant derivative in 3-D space. 

• cfJ : 3-D state of stress ( i j=l, 2, 3). 

•P 1 : external volume force. 

• u a : Gaussian coordinates of surface ( ct=l, 2). 

• z : normal coordinate of surface. 

and applying the following Kirchhoff-Love hypothesis to virtual displacements 

fo a (u°,z) -(«' -zd[) *(6v r - zdq y ) 
dv 3 (u a ,z ) - 6w 

where 


(6.2) 

(6.3) 


a r 

• « : component of mixed metric tensor. 

d r 

* a : component of mixed curvature tensor. 

Here the rotation dq y is defined as 

(64) 

Then, following (he lengthy derivations by F. I. Niordsonin his Shell Theory [I], the principle 
of virtual work gives 

dK'j ]dA 

~ S5[F a 6v a +pdw]dA + $jr a 6v a +M°fr a + Qdw ]ds (6.5) 

Where 
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• N a P : effective membrane stress tensor. 

• : effective moment tensor. 

• dE^p : virtual strain tensor. 

• dK a p : virtual bending tensor. 

• F 0 , p . effective load components in surface and normal directions, respectively. 

• T 01 •' membrane force vector acting on the boundary. 

M . moment vector acting on the boundary. 

• Q . supplemented shear force on the boundary. 

• dv a , dw : virtual displacements in plane and transverse directions, respectively. 

• dr a =e a Pdqp 

• e a ^ : alternate tensor. 


Mathematically, the strain tensor E a p and bending tensor are defined as 

E ap = l( a 'a0 - a ae ) ( 6 . 6 ) 

K = d" - d _ 

°P op afi (6.7) 


where the superscript asterisks indicate the deformed state and can be derived as follows : 

a 'ofS = % + +/V + P°P 0r + q ° q e ( 

1 

d 'a, - (fr)’t(l + ^ +p/a)(d t0 +D p q a + d\ P<Tf ) 

~{q p ♦e"** q r p^)(D pPap -d 0p q a )] (( 


The generalized two-dimensional displacement gradient P and its determinant in 
equation (6.9) are expressed as : 


Pa,, ~ 

V - det (fV 


( 6 . 10 ) 

( 6 . 11 ) 
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After linearization, the strain and bending tensors can be expressed in terras of 
displacements as follows : 


+ D 0 V a ) - d^W 


~Dn 


afl 


+ pV r + d fif D a v r + v r - d ^d r a w 


( 6 . 12 ) 
(6.13) 


be 


Considering the isotropic thin elastic shell for simplicity. The 


N °* v)E* + va* e\] 


M 




Eh 




-t(l - v )Ar^ + va aP K r r ] 


constitutive equations will 

(6. 14) 

(6.15) 


The substitution of the last four equations into equation (6.5) will result in the weak 
forra whtch ts expressed in terras of dtsplaceraent vectors. Before applying the finite eleraen, 
method to solve equatron (65), a univeisal methodology is outlined ra the next paragraph. 

6.3 Methodology for solving shell problem by FEM 

Afttr ,he Potions of mathematical formulation, i, ,s necessary to discretize equatton 
(65) to solve shell problems by FEM. However, as equations (612) to (615) show the 
calculations of const, tutive equations and stram-displacemen. telations mvolve the evaluates 
erf raetnc, curvature ransom, and covenant derivatives. Moreover, in general, the calculations 
o covanant denvattves require the computations of Christoffel symbols. If a given geometric 
omain is complex, these calculations will be tedious. With the help of symbolic and algebraic 
manipulation, these tough tasks can be performed by simply giving ihe parametric equations of 
the surface. Here the outline of methodology to solve shell problems is presented as follows : 

1 Finding the parametric equations of the middle surface for a given shell. 

2 Calculating the Christoffel symbols, the metric and curvature tensors based on the 
paramedic equations. If the parametric equations are chosen correctly, the resulting 
metnc and curvature tensor should obey the integrabilily condition. In other words, they 

should not violate the Coddazzi-Mainardi equations. Gaussian equations and the regular 
condition. ° 


3. SubsUtuting the metric, curvature tensors, and Christoffel symbols into constitutive 
equations and strain-displacement equations. 
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4 Discretizing the variational form of equation (6.5) and solving it by finite element 
method. 


The above methodology is universal for any shell problems. Different shell geometries 
can be solved in the same way by simply feeding appropriate parametric equations into the 
symbolic and algebraic manipulator REDUCE. Based on this methodology, an example of 
plate problem is shown in the following paragraph. 

6.4 Symbolic and algebraic manipulation application to plate problems 
6.4.1 Methodology 

1 Three parametric equations are chosen as follows : 
f 1= u ! , f 2 =u 2 , f 3 =constant 

2 Calculate surface metric, curvature tensors, and Christoff el symbols symbolically. These 
calculations are based on their fundamental definitions, given by: 

• metric tensor : 


i /» ^ 

^ afi ~ f af 0 


(6.16) 


• curvature tensor : 


(6.17) 

where / , , / 2 and X* are surface Gaussian and normal coordinates. They are shown in 
Figure 5.2 pictorially, and X* is defined by : 


~7 W r k 

a e JJ.i 


(6.18) 


• The 2nd kind of Christoffel symbols is defined as: 




(6.19) 


The REDUCE program for calculating equations (6.16), (6.17), (6.18) and (6.19) is 
presented as follows. The results follow the program. 
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% ******** * ** * ** * * * ****** * ** * *#** * * * * * * * * * * * * * * * * * * * * * * * * $ * * £ £ 
% REDUCE PROGRAM FOR CALCULATING METRIC, CURVATURE 
% TENSORS & CHRISTOFFEL SYMBOLS 

%*********************************3Kt**********t*********^*** 


% - 

%INPUTTING THE PARAMETRIC EQUATIONS 


MATRIX A(2,2),CA(2,2),F(2,3),D(2,2); 

ARRAY X(3),C1(2,2,2),C2(2,2,2),N(3),E(3,3,3),U(2); 

U(2):=P; 

OFF PERIOD; 

X(I):=U(1); 

X(2):=U(2); 

X(3):=CONSTANT ; 

FOR I:=l;2 DO FOR J:=i:3 DO 
F(I,J):=DF(X(J),U(I)); 

FOR ALL T 1 LET COS(T 1 )**2+SIN(T 1)**2= 1 ; 

% - 

%CALCULATING COVARIANT METRIC TENSOR 

% — 

FOR M:=l:2 DO FOR N:=l:2 DO 

A(M,N):=FOR I:=l:3 SUM DF(X(I),U(M))*DF(X(I),U(N)); 

A* = A* 

DETA:=DET(A); 


% 

%C A LCULATING CONTRA VARIANT METRIC TENSOR 

% 

FOR L:=l:2 DO FOR M:=l:2 DO 
IF L=1 AND M=1 THEN CA(L,M):=A(2,2)/DETA 
ELSE IF L NEQ M THEN CA(L,M):=-A(M,L)/DETA 
ELSE CA(L,M):=A( 1,1 )/DETA; 


%~ - 

%CALCULATING THE 1ST & 2ND CHRISTOFFEL SYMBOL 

% 

WRITE "THE 2ND CHRISTOFFEL SYMBOL" • 

FOR L:=l:2 DO FOR M:=l:2 DO FOR N:=l:2 DO 
C 1(L,M,N):=( 1/2)*(DF(A(L,N),U(M))+DF(A(M,N),U(L)) 
-DF(A(L,M),U(N))); 

FOR L:=l:2 DO FOR M:=l:2 DO FOR N:=l:2 DO 
«C2(L,M,N):=FOR I:=l:2 SUM CA(L,I)*C1(M,N,IV 
WRITE "C2(",L,",",M,",",N,")=",C2(L,M,N)»; 


% - 

%CALCULATING ALTERNATING TENSOR E(U,K) 

% 

FOR I:=l:3 DO FOR J:=l;3 DO FOR K:=l:3 DO 
«IF (I=J OR J=K OR I=K) THEN E(I,J,K):=0 
ELSE IF (1=1 AND J=2 AND K=3) OR (1=2 AND J=3 AND K=l) OR 
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(1=3 AND J=1 AND K=2) THEN E(1,J,K):=1 
ELSE E(I,J,K):=-1; WRITE E(I,J,K)»; 

%- - 

%CALCULATING NORMAL VECTOR 

%- - - 

FOR I:=l:3 DO 

N(I):=(FOR J:=l:3 SUM FOR K:=l:3 SUM 
E(I,J,K)*F( U)*F(2,K))/SQRT(DETA); 

% 

%CALCULATING CURVATURE TENSOR 
----- - - — 

FOR J:=l:2 DO FOR K:=l:2 DO 
D(J,K):=FOR I:=l:3 SUM N(I)*DF(F(J,I),U(K)); 
FOR ALLT1 CLEAR COS(Tl)**2+SIN(Tl)**2; 
FOR ALL T1 LET COS(Tl)**2-l=-SIN(Tl)**2; 


% - 

%OUTPUTTING RESULTS 

% 

OFF PERIOD; 

OFF ECHO; 

OUT "SURFACE"; 

WRITE " 

WRITE " THE COMPONENTS OF METRIC TENSOR " ;’ 

WRITE " - — "• 

A:=A; 

WRITE " ----- » 

WRITE " THE COMPONENTS OF CURVATURE TENSOR" 

WRITE " 

D:=D; 

WRITE " - - "; 

WRITE " THE CHRISTOFFEL SYMBOLS"; 

WRITE " - - — 

FOR I:=l:2 DO FOR J:= 1:2 DO FOR K:=l:2 DO 
WRITE "CHRIS(",I,",",J,",",K,")=",C2(I,J,K); 

SHUT "SURFACE"; 

BYE; 


The outputs of REDUCE are presented as follows : 


THE COMPONENTS OF METRIC TENSOR 


A(l,l) := 1 
A(l,2) := 0 
A(2,i) := 0 
A(2,2) := 1 


THE COMPONENTS OF CURVATURE TENSOR 


D(l,l) := 0 
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D(l,2) := 0 
D(2,l) := 0 
D(2,2) := 0 


THE CHRJSTOFFEL SYMBOLS 


CHRIS(1,1,1)=0 
CHRJS(1,1,2)=0 
CHRIS(1,2, 1)=0 
CHRIS( 1,2,2)=0 
CHRIS(2,1,1)=0 
CHRIS(2,1,2)=0 
CHRIS(2,2, 1)=0 
CHRIS(2,2,2)=0 


As the output from REDUCE shows, all of the Chnstoffel symbols vanish in the plate 

case. According to the above solutions, the metric and curvature tensors in matrix form 
are 

1 ■ (o *5 • K# 1 - ( 0 (6.20) 

3. After substituting the calculated metric and curvature tensors, the constitutive equations 
are as follows 


• For membrane : 


V 

N 12 

Eh 

1-v * 

-1 0 v 

0 1 - v 0 


E 11 ■ 

iV 22 


y 0 l 


E u 


•For bending : 


'M n ' 

1 

■i 

0 

V 


k" 

M n 

Eh 

0 

1 - V 

0 



laci-v 1 ) 


K 

M 12 


V 

0 

1 


K ” 


( 6 . 21 ) 


( 6 . 22 ) 


The relationships between linearized strain-displacement and bending curvature tensors 
are from equations (6. 12) and (6. 13) : 


E c 4 > “ k(D a v $ +D p v a ) - a + v a ^ ) 
K - D X D w - w 

.afi 


(6.23) 

(6.24) 
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where the disappearance of Chnstoffel symbols in this case eliminates the distinction 
between covariant different, ation and ordinary partial differentiation. In add.t.on, the 

umty of the metric tensors remove the distinction between contravariant and covariant 
tensors. 


4. Finite element method starts (see the following paragraph). 

6.4.2 Finite element discretization 

, rsT , ^ e ‘ ement Ch ° Sen for ,his *°P ic ' s combination of the constant strain Wangle 
(CST) element with the Cheung, King, Zientnewicz (CKZ) triangle element The 

considerations of selecting this element will be discussed later. The shape functions for this 
element are 

L a “I c 

N<x ” £cr + (to “ ~ % r ) 

N '“ ’ 2A[C ' ( 5”5, + T?,5,5,) - c „( + £,?,)) (6 25) 

-2AI + 11,5,1,) - ft, (5^, + 

Where (a,p,y) is the permutation of (1,2,3) and no summation convention is applied 
equation (6.25). The constants b ( c-, and A are defined as follows : 


b = 

h - Ll Zi 

1 2 A * 

^2 2 A 

x l“ x i 

r =as — 


C l 2 A » 

C Z ** 2A ’ 

2A = _r 2 y 3 - 

*3?2 + * 3 ? 


V i - y 
^ 3 = " 2 A 


2A 


(6.26) 


Then the deflections u, v, w in local coordinates jt, y, z direction can be interpolated by 

3 

l -l 

The substitution of equations (6.21) to (6.27) into (6.5) gives the discretization form 

-ff/dA + fs r ds + fM r ds ( 6 . 28 ) 


(6.27) 
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Where 




D . 


Eh 

I2(l-v ) 


1 0 V 

0 if 0 

v 0 1 

j 

v 
0 

0 1 


ri 


o 

0 if 


(6.29) 


(6.30) 


B. 


c x 

0 


0 0 


0 

b. 


0 
b . 


0 
l » 


o 0 0 0 L. 

0 L, 0 0 0 0 


r 0 
r 0 
r 0 


0 


j 

■5. 


0 — 

0 — 

0 0 0 0 I 3 

^,0000 


0 0 0 0 
£,000 


(6.31) 


'*i b 2 b 3 0 0 0 

fl 2 “ C 1 c 2 C 3 *i *2 *3 

0 0 0 c. c c 


r * 
*■ 
a 

*7 

<» 

*7 


o 

0 

0 


0 — 
0 — 

0 — 

■5. 


*>, b. 




(6.32) 


M - [A/ 1 A/ 2 ] 


c > C 2 C 3 

b t b j 


_£_i 

<? 

«7 

iT 


AT 


(6.33) 
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N. 
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0 
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6.4.3 Numerical results and post-process 


Three boundary conditions are applied to 


N - N >, 0 0 N , N„] 

“V, ”, >%, »,,,l 

( 6 . 34 ) 

( 6 . 35 ) 

( 6 . 36 ) 

0 0 L 3 0 0 o 0 

0 0 o l 3 0 0 0 

0 0 N, N } ' N}) 

( 6 . 37 ) 


test problem. They are : 


1. In-plane umax.al tension [see Figure 6.1] ... In this case, only the membrane component 
contribute to the stiffness matrix. The consistent load vector is due to the boundary 
traction S only. The REDUCE programs to the stiffness and the consistent load vector 

are presented as follows. In addition, the stress distribution can also be calculated 
symbolically. 
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%* ****************** ********************* ***##**##**#*****## 

% Symbolical program for making stiffness matrix for plate tension problem 
%**************************** ****************** *#*####*##**# 
matrix bc(3,6),In(6,15),d(3,3),b(3,15),s(15,15); 
array n(3); 

n( l):=sl ;n(2):=s2;n(3):=s3; 

In(l,l):=l;ln(2,6):=l;ln(3,ll):=l;ln(4,2):=l; 
ln(5,7):= 1 ;ln(6, 12):= 1 ; 

bc:=mat((bl,b2,b3,0,0,0),(cl,c2,c3,bl,b2,b3),(0,0,0,cl,c2,c3)); 

d:=mat(( 1,0, v), (0,(1 -v)/2,0),(v,0,l)); 

b:=bc*ln; 

s:=tp(b)*d*b*pj*ye*h/(2*(l-v**2)); 

b 1 : =( y 2- y3 ) /pj ; b2: =( y3 - y 1 )/pj ;b3 : =( y 1 -y2)/pj ; 

c 1 : =( x3 -x2)/pj ;c2 : =(x 1 -x3 )/pj ;c3 : =(x2- x 1 )/pj ; 

for i:=l: 15 do «a:=for j:=l: 15 sum s(i,j);write a»; 

for i:=l: 15 do «c:=for j:=l: 15 sum s(j,i);write c»; 

on fort; 

off period; 

off echo; 

out "plane, ftn"; 

write " subroutine ske(xl,yl,x2,y2,x3,y3,s) H ; 

write" dimension s(15,15)"; 

write " implicit real*8(a-h,o-z)"; 

write " Pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3)"; 

write" ye=2.1*1.0e07"; 

write " v=0.29"; 

write " h=0.2"; 

for i:=l: 15 do for j:=l: 15 do 

if j>=i then write " s(\i,",",j,")=",s(i,j) 

else write * s(",i,",",j,")=s(",j,",",i,") H ; 

write " return"; 

write" end"; 

shut "plane. ftn"; 

bye; 

% Symbolic program to calculate the load vector for plate tension problem. 

(fo********************************************************* 

matrix f(l,2),sn(2,15),fn(15,l); 
array n(3),ff(15); 
n(l):=sl;n(2):=s2;n(3):=s3; 
f:=mat((fl,f2)); 

for i:=l step 5 until 11 do«sn(l,i):=n((i-l)/5+l); 

sn(2,i+l):=n((i-l)/5+l)»; 

fn:=tp(f*sn); 

s2:=0; 

s3:=l-sl; 

for i:= 1: 15 do «al:=int(fn(i,l),sl);a2:=sub(sl=l,al)-sub(sl=0,al); 

fn(i,l):=a2*rl*h»; 

for i:=l;15do ff(i):=fn(i,l); 

on fort; 

off echo; 

off period; 
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out "15load.ftn"; 

write " subroutine load(xl,yl,x2,y2,x3,y3,fl,f2,f3,fe)"; 

write " implicit real*8(a-h,o-z)"; 

write " dimension fe( 15)"; 

write " rl=sqrt((xl-x3)**2+(yl-y3)**2)"; 

write " h=0.2"; 

for i: = l: 15 do write " fe(",i,")=",ff(i); 

write " return"; 

write" end"; 

shut "151oad.ftn"; 

bye; 

%**** ********************* ****************** 

% REDUCE program to construct subroutine to compute 
% stress distribution for plate tension problem. 
%********************** ********* ************ 
matrix d(3,3),bc(3,6),fn(6,l),stre(3,l); 
array n(3) ; 

n(l):=sl;n(2):=s2;n(3):=s3; 

bc:=mat((bl,b2,b3,0,0,0),(cl,c2,c3,bl,b2,b3), (0,0,0, cl,c2,c3)); 
d:=mat(( l,0,v),(0,( l-v)/2,0),(v,0, 1)); 
operator u; 

for i: = l:6 do fn(i,l):=u(i); 
stre:=d*bc*fn*ye/( l-v**2); 
bl:=(y2-y3)/pj;b2:=(y3-yl)/pj;b3:=(yl-y2)/pj; 
c I:=(x3-x2)/pj ;c2:=(x 1 -x3)/pj ;c3:=(x2-x 1 )/pj ; 
on fort; 
off echo; 
off period; 
out "st.ftn"; 

write " subroutine stress(xl,yl,x2,y2,x3,y3,u,st)"; 
write" implicit real*8(a-h,o-z)"; 

write " dimension u(6),st(3)"; 

write " v=0.29"; 

write " pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3)"; 

write" ye=2. I*1.0e07"; 

for i:=l:3 do write " st(",i,")=",stre(i,l); 

write " return"; 

write" end"; 

shut "st.ftn"; 

bye; 


The resultant fortran subroutines obtained from the above programs to compute stiffness 
matrix, load vector and stress distribution are presented as follows. Since the expression 
of stiffness matrix is quite lengthy, only a part of it is showed. The interested 
researchers may refer to the PH. D. thesis of Wen-LangTsai [51] for details. 

SUBROUTINE SKE(X1,Y 1,X2,Y2,X3,Y3,S) 

IMPLICIT REAL*8(a-h,oz) 

DIMENSION S( 15, 15) 

Pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3) 

YE=2. I* 1.0E07 
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V=0.29 

H=0.2 

s(l,l)=(H*YE*(V*X2**2-2*V*X2*X3+V*X3**2-X2**2+2*X2* 

. X3-X3**2-2*Y2**2+4*Y2*Y3-2*Y3**2))/(4*PJ*(V**2-1)) 
s(l,2)=(H*YE*(V*X2*Y2-V*X2*Y3-V*X3*Y2+V*X3*Y3+X2*Y2- 
. X2*Y3-X3* Y2+X3*Y3))/(4*PJ*( V**2- 1)) 
s(l,3)=0 
s(l,4)=0 
s( 1 ,5)=0 

s(l,6)=-(H*YE*(V*Xl*X2-V*Xl*X3-V*X2*X3+V*X3**2-Xl*X2 
. +X1*X3+X2*X3-X3**2+2*Y2*Y3-2*Y2*Y1-2*Y3**2+2*Y3*Y1) 

. )/(4*PJ*(V**2-l)) 

s(l,7)=-(H*YE*(2*V*Xl*Y2-2*V*Xl*Y3+V*X2*Y3-V*X2*Yl-2 
* V*X3*Y2+V*X3*Y3+V*X3*Y 1-X2*Y3+X2*Y 1+X3*Y3-X3*Y 1)) 
. /(4*PJ*(V**2-1)) 
s(l,8)=0 
s( 1 ,9)=0 
s(l,10)=0 

s(l,l 1)=(H*YE*(V*X1*X2-V*X1*X3-V*X2**2+V*X2*X3-X1*X2 
. +X 1 *X3+X2**2-X2*X3+2*Y2**2-2*Y2*Y3-2*Y2*Y 1+2*Y3*Y 1) 

. )/(4*PJ*(V**2-l)) 

s(l,12)=(H*YE*(2*V*Xl*Y2-2*V*Xl*Y3-V*X2*Y2+2*V*X2*Y3 
-V*X2*Y 1- V*X3*Y2+V*X3*Y 1-X2*Y2+X2*Y 1+X3*Y2-X3*Y 1)) 

. /(4*PJ*(V**2-1)) 
s(l,13)=0 
s(l,14)=0 
s(l,15)=0 
s(2,l)=s( 1,2) 

s(2,2)=(H*YE*(V*Y2**2-2*V*Y2*Y3+V*Y3**2-2*X2**2+4*X2 
. *X3-2*X3**2-Y2**2+2*Y2*Y3-Y3**2))/(4*PJ*(V**2-1)) 
s(2,3)=0 
s(2,4)=0 
s(2,5)=0 

s(2,6)=(H*YE*(V*Xl*Y2-V*Xl*Y3+2*V*X2*Y3-2*V*X2*Yl-V* 

. X3*Y2-V*X3*Y3+2*V*X3*Y1-X1*Y2+X1*Y3+X3*Y2-X3*Y3)) 
/(4#pj*(V**2-l)) 

s(2,7)=(H*YE*(V*Y2*Y3-V*Y2*Y1-V*Y3**2+V*Y3*Y1+2*X1* 

. X2-2*X1*X3-2*X2*X3+2*X3**2-Y2*Y3+Y2*Y 1+Y3**2-Y3*Y 1)) 

. /(4*PJ*(V**2-1)) 
s(2,8)=0 
s(2,9)=0 
s(2,10)=0 

s(2, 1 1)=-(H*YE*(V*X 1 * Y2- V*X 1 * Y3+V*X2* Y2+ V*X2* Y3-2* V* 

. X2*Y 1-2*V*X3*Y2+2*V*X3*Y 1-X1*Y2+X1*Y3+X2*Y2-X2*Y3)) 
. /(4*PJ*(V**2-1)) 

s(2,12)=-(H*YE*(V*Y2**2-V*Y2*Y3-V*Y2*Yl+V*Y3*Yl+2*Xl 
. *X2-2*X1*X3-2*X2**2+2*X2*X3-Y2**2+Y2*Y3+Y2*Y 1-Y3*Y1 
. ))/(4*PJ*(V**2-l)) 


s(15, 15)=0 

return 

end 
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subroutine load(x 1 ,y 1 ,x2,y2,x3,y3,f 1 ,f2,f e) 

IMPLICIT REAL*8(a-h,o-z) 

dimension fe(15) 

rl=sqrt((x3-x2)**2+(y3-y2)**2) 

h=0.2 

fe( 1)=0 

fe(2)=0 

fe(3)=0 

fe(4)=0 

fe(5)=0 

fe(6)=(H*Fl *RL)/2 
fe(7)=( H* F2* RL)/2 
fe(8)=0 
fe(9)=0 
fe(10)=0 

fe( 1 l)=(H*Fl*RL)/2 

fe( 12)=(H*F2*RL)/2 

fe(13)=0 

fe(14)=0 

fe(15)=0 

return 

end 

subroutine stress(x 1 ,y 1 ,x2,y2,x3,y3,u,st) 
implicit reaJ*8(a-h,o-z) 
dimension u(6),st(3) 
v=0.29 

pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3) 

YE=2.1*1.0E07 

st(l)=(YE*(U(6)*V*Xl-U(6)*V*X2-U(5)*V*Xl+U(5)*V*X3+U 
(4)*V*X2-U(4)*V*X3+U(3)*Y2-U(3)*Y 1-U(2)*Y3+U(2)*Y 1- 
. U(1)*Y2+U(1)*Y3))/(PJ*(V**2-1)) 
st(2)=-(YE*(U(6)*V*Y2-U(6)*V*Yl-U(6)*Y2+U(6)*Yl-U(5) 

. * V*Y3+U(5)* V*Y 1+U(5)*Y3-U(5)*Y 1-U(4)*V*Y2+U(4)* V*Y3 
. +U(4)*Y2-U(4)*Y3+U(3)*V*X1-U(3)*V*X2-U(3)*X1+U(3)* 

. X2-U(2)*V*X1+U(2)*V*X3+U(2)*X1-U(2)*X3+U(1)*V*X2-U( 
1)* V*X3-U( I)*X2+U( 1)*X3))/(2*PJ*(V**2- 1)) 
st(3)=(YE*(U(6)*Xl-U(6)*X2-U(5)*Xl+U(5)*X3+U(4)*X2-U 
. (4)*X3+U(3)*V*Y2-U(3)*V*Y1-U(2)*V*Y3+U(2)*V*Y1-U(1) 

. *V*Y2+U(1)*V*Y3))/(PJ*(V**2-1)) 
return 
end 


The stress distribution of plate under tension is plotted by PATRAN in Figure 6.2. The 
stress concentration is visible in the top of the hole. One interesting phenomenon that 
should be mentioned here is that the stress at the middle top of the plate is less than the 
applied load. This is contributed by the bending effect which produces the compression 
in the top fiber of the plate. 
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no 


Figure 6.2 : Stress distribution of plate under uniaxial tension 


2. Four-edge simply supported plate bending with uniform transverse load [see Figure 
6.3] — In this case, the membrane component of stiffness matrix is neglected. The 

REDUCE program to make the bending component of stiffness matrix is shown as 
follows : 

% *************^ 

ng P late ^'ng stiffness matrix 

array n(9); 

"T)^= s f "ft ^ bc < 3 - 6 >-*—<3^) -<«3.3 ) ,ske(9,9),s(9,9); 

n(2):=pj*(c3*(s 1 * *2*s2+s 1 *s2*s3/2)-c2*(sl **2*s3+sl *s2*s3/2» • 
n(3):-pj*(b2*(sL**2*s3+sl*s2*s3/2)-b3*(sl**2*s2+sl*s2*s3P)V 
n(4):=s2+s2*s3*(s2-s3)+s2* s i*(s2-sl); }) ’ 

n(5):=pj*(cl*(s2**2*s3+sl*s2*s3/2)-c3*(s2**2*sl+sl*s2*s3/2)V 
n(6):=pj*(b3*(s2**2*sl+sl*s2*s3/2)-bl*(s2**2*s3+sl*s2*s3/2Y>'- 
n(7):=s3+s3* s i*(s3-sl)+s3*s2*(s3-s2)- }) ’ 

n(8):= PJ *(c2*(s3**2*sl+sl*s2*s3/2)-ci*(s3**2*s2+sl*s2*s3/2))- 
n(9):=pj*(bl*(s3**2*s2+sl*s2*s3/2)-b2*(s3**2*sl+sl*s2*s3/ :: >l') , • 
for i:=l:9do«sn(l,i):=df(n(i),sl);sn(2,i):=df(n(i),s2)- " ’ 

sn(3,i):=df(n(i),s3)»; 
bc:=mat((bl,b2,b3),(cl,c2,c3)); 
bn:=bc*sn; 

f or i: = l : 9 do <<ssn(l,i) : =d f ( bn (l,i) > s l);ssn(2,i):=df(bn(l,i),s2); 
ssn(3,i):=df(bn(l,i),s3);ssn(4,i):=df(bn(2,i),sl); 
ssn(5,i):=df(bn(2,i),s2);ssn(6,i):=df(bn(2,i),s3)>>- 

s3:=l-sl-s2; 

te^=lx^*ssn b2,b3,0,0,0),(Cl ’ C ^ ,C ^ ^’b3)’(0’0.0,cl,c2,c3)); 

D: =M AT ((1 ,0, V) ,(0,( 1 - V)/2,0) ,( V,0, 1 )) ; 

ske:=tp(bssn)*d*bssn*pj$ 

for i:=l:9 do for j:=l:9 do 

if j>=i then «tem:=int(ske(ij),s2); 

tern 1 :=sub(s2= 1 -s 1 ,tem)-sub(s2=0,tem) ; 

tem3:=int(teml,sl); 

s(i,j):=(sub(sl=l,tem3)-sub(sl=0,tem3))*ye*h**3/(l'>*(l- v**2)I» 

else s(i,j)=s(j,i); 

b 1 : =(y2-y3)/pj ;b2: =(y3-y 1 )/pj ;b3:=(y 1 -y2)/pj ; 

c 1 : =( x3 -x2)/pj ;c2:=(x 1 -x3)/pj ;c3 : =(x 2 - x 1 )/pj ; 

on fort; 
off echo; 
off period; 

CARDNO!*:=10; 
out "z.ftn"; 

WRITE " SUBROUTINE SKE1 (XI Y1 X2 Y2X3 Y3 SV- 
WRITE " IMPLICIT REAL*8(A-H O-zf ’ ’ } ’ 

WRITE” DIMENSION S(2, 9)"; 

WRITE " Pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3)"; 

WRITE” YE=2. 1* 1.0E07”; 

WRITE" V=0.29”; 

WRITE" H=0.2"; * 


FOR I:=l:2 DO FOR J:=l:9 DO 
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IF j>=i THEN WRITE " S(M,",",J,")=",S(I,J) 

ELSE WRITE " S(",I,",",J,")=S(",J,",",I,T; 

WRITE" RETURN"; 

WRITE" END"; 

WRITE " SUBROUTINE SKE2(X1,Y 1,X2,Y2,X3,Y3,S)"; 
WRITE " IMPLICIT REAL*8(A-H,0-Z)"; 

WRITE " DIMENSION S(4,9),sl(2,9)"; 

WRITE " pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3)"; 

WRITE" YE=2.1*1.0E07"; 

WRITE" V=0.29"; 

WRITE" H=0.2"; 

WRITE " CALL SKE1(X1,Y 1,X2,Y2,X3,Y3,S1)"; 
WRITE" DO 20 1=1,2"; 

WRITE" DO 20 J= 1,9"; 

WRITE" 20S(I,J)=S1(U)"; 

FOR I:=3:4 DO FOR J:=l:9 DO 

IF j>=i THEN WRITE " S(M,",\J,")=",S(I,J) 

ELSE WRITE " S(",I,",",J,")=S(",J,",",L")"; 

WRITE" RETURN"; 

WRITE" END"; 

WRITE * SUBROUTINE SKE3(X1,Y 1,X2,Y2,X3,Y3,S) M ; 
WRITE ” IMPLICIT REAL*8(A-H,0-Z)"; 

WRITE " DIMENSION S(6,9),s2(4,9)"; 

WRITE " pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3)"; 

WRITE" YE=2. 1* 1.0E07"; 

WRITE" V=0.29"; 

WRITE" H=0.2"; 

WRITE " CALL SKE2(X1,Y1,X2,Y2,X3,Y3,S2)"; 
WRITE" DO 20 1=1,4"; 

WRITE" DO 20 J= 1,9"; 

WRITE" 20S(U)=S2(I,J) n ; 

FOR I:=5:6 DO FOR J:=l:9 DO 

IF j>=i THEN WRITE " S(",I,7\J,")=",S(U) 

ELSE WRITE " S(",I,",",J,")=S(",J,",",I,")"; 

WRITE" RETURN"; 

WRITE" END"; 

WRITE " SUBROUTINE SKE(XI,Y1,X2,Y2,X3,Y3,S)"; 
WRITE " IMPLICIT REAL*8(A-H,0-Z)"; 

WRITE " DIMENSION S(9,9),s3(6,9)"; 

WRITE " pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3)"; 

WRITE" YE=2. 1*1.0E07"; 

WRITE" V=0.29"; 

WRITE" H=0.2"; 

WRITE " CALL SKE3(X1,Y 1,X2 > Y2,X3,Y3,S3)"; 
WRITE" DO 20 1=1,6"; 

WRITE" DO 20 J= 1,9"; 

WRITE" 20 S(I,J)=S3(I,J)"; 

FOR I:=7:9 DO FOR J:=l:9 DO 

IF j>=i THEN WRITE " S(",I,",",J,")=",S(I,J) 

ELSE WRITE " S(",I,",",J,")=S(",J,",",I,T; 

WRITE" RETURN"; 

WRITE" END"; 

SHUT "z.ftn"; 
bye; 
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The resultant fortran subroutine is too large (around 135 pages) to be presented 
completely here. The following is only a small portion of it. 


SUBROUTINE SKE1(X1,Y1,X2,Y2,X3,Y3 S) 
IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION S(l,9) 
pj=(xl-x3)*(y2-y3)-(yl-y3)*(x2-x3) 

YE=2. 1* 1.0E07 

V=0.29 

H=0.2 


ANS5=-4*Y3*Y 1**3+2*Y 1**4 

ANS4=-4*X2**2*Y2*Y 1+8*X2**2*Y3**2+2*X2**2*Y1**2-16* 
X2*X3**3-16*X2*X3*Y2**2+32*X2*X3*Y2*Y3-16*X2*X3*Y3 
. **2+5*X3**4+8*X3**2*Y2**2-16*X3**2*Y2*Y3+10*X3** , ' > * 

. Y3**2-4*X3**2*Y3*Y1+2*X3**2*Y1**2+5*Y2**4-16*Y2**3* 

. Y3-4*Y2**3*Y1+24*Y2**2*Y3**2+6*Y2**2*Y1**2-16*Y*"' , *Y3 
• **3-4*Y2*Y 1**3+5*Y3**4-4*Y3**3*Y 1+6*Y3**2*Y 1 **?+ 

. ANS5 

ANS3=2*Xl**4-4*Xl**3*X2-4*Xl**3*X3+6*Xl**2*X2**2+6* 
X1**2*X3**2+2*X1**2*Y2**2-4*X1**2*Y2*Y1+2*X1**" > *Y3 
**2-4*Xl**2*Y3*Yl+4*Xl**2*Yl**2-4*Xl*X2**3-4*Xl*X' ) * 
Y2**2+8*X1*X2*Y2*Y 1-4*X1*X2*Y l**2-4*Xl*X3**3-4*xI"* 

. X3*Y3**2+8*X1*X3*Y3*Y 1-4*X1*X3*Y1**2+5*X2**4-16*X'> 

. **3*X3+24*X2**2*X3**2+10*X2**2*Y2**' , -16*X2** , ‘ , *Y'>*Y3 
. +ANS4 “ “ 


ANS2=H**3*YE*ANS3 
ANS1=ANS2/(18*PJ**3*(V**2-1)) 
S(l,l)= -ANSI 
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The load vector in this case is for a uniform transverse pressure only. The REDUCE 
program and its output are shown as follows : 


%*******************%**%*%*%%%*%%%%% 

% REDUCE program to construct the load vector 
% for plate bending problem. 

array n(9),fe(9); 

n(l):=sl+sl*s2*(sl-s2)+sl*s3*(sl-s3); 

n(2):=pj*(c3*(sl**2*s2+sl*s2*s3/2)-c2*(sl**2*s3+sl*s2*s3/2)V 

n (3):=pj*(b2*(sl**2*s3+sl*s2*s3/2)-b3*(sl**2*s2+sl*s2*s3/2)V 

n(4):=s2+s2*s3*(s2-s3)+s2*sl*(s2-sl); 

n(5):=pj*(cl*(s2**2*s3+sl*s2*s3/2)-c3*(s2**2*sl+sl*s2*s3/2))- 

n(6):=pj*(b3*(s2**2*sl+sl*s2*s3/2)-bl*(s2**2*s3+sl*s2*s3/2)V 

n(7):=s3+s3*sl*(s3-sl)+s3*s2*(s3-s2); 

n(8):=pj*(c2*(s3**2*sl+sl*s2*s3/2)-cl*(s3**2*s2+sl*s2*s3/2)V 

n(9):=pj*(bl*(s3**2*s2+sl*s2*s3/2)-b2*(s3**2*sl+sl*s2*s3/2)V 

s3:=l-sl-s2; 

for i:=l:9 do «teml:=int(n(i),s2); 

tem2:=sub(s2= 1 -s 1 .tern 1 )-sub(s2=0,tem 1); 
tem3:=int(tem2,sl); 

fe(i):=(sub(sl=l,tem3)-sub(sl=0,tem3))*pj*f3; 
write i,fe(i)»; 

b 1 : =(y2-y3)/pj ;b2: =(y3-y 1 )/pj ;b3: =( y 1 -y2)/pj ; 
Cl:=(x3-x2)/pj;c2:=(xl-x3)/pj;c3:=(x2-xl)/pj; 

off period; 
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off echo; 

ON FORT; 
out "zload.ftn"; 

write " SUBROUTINE LOAD(Xl Yl X2 Y2 Yi yi P 3 trc\i* 
write * IMPLICIT REAL*8(A-H aZ)-’ ’ Y2,X3,Y3,R,FE) : 

WRITE ■ DIMENSION FE(9)" ' ’ 

" Pj=(xi-x3)*(y2-y3)-(yl-y3)*(x2-x3i"- 

FOR I:= 1 :9 DO WRITE " FET I FFnl ’ 

WRITE " RETURN"- * ” ,FE(I); 

write " end"; 
shut "zload.ftn"; 
bye; 

SUBROUTINE LOAD(Xl,Yl X2 Y? X3 Yl Fl pp^ 
IMPLICIT REAL*8(A-H O-Z) ’ ’ ,FE) 

DIMENSION FE(9) ’ 

lwg 3M,, - ,3w ' ,3) 

FE(2)=-(PJ*F3*(2*X l-X2-X3))/48 
FE(3)=(PJ*F3*(Y2+Y3-2*Yl))/48 
FE(4)=(PJ*F3)/6 )} 

FE(5)=(PJ*F3*(Xl-2*X2+X3))/48 
FE(6)=-( PJ* F3 * (2* Y 2- Y3- Y 1 ) )/48 
FE(7MPJ*F3)/6 ” 

FE(8)=(PJ*F3*(Xl+X2-2*X3))/48 
FE(9)=(PJ*F3*(Y2-2*Y3 + Y1))A« 
return }) 

end 


The deformed shape ,s shown i„ Figure 6.4 and ihe stress disinbution pattern ,s 
Ftgure 6.5. In addttion. three dtfferen, stzes of mesh are tested to tnvesttga 
convergence of solutton. They are shown in Figure 6.6. The convergence tend * 
presented m Ftgure 6.7 which is the plot of error vs, element mesh size. 

3- Four-edge clamped plate bending under un.foon load - the only dtfference between this 
ase an t e last case is the boundary constraints. The extra slope constraints are 

“ " t “ - -Pected dial ihe solutton wil, be sZ 

hm, that of the stmply-supported case. The deformed shape is shown in Ftgure 6 8 The 
iffer phenomenon is visible by comparing Figures 6.4 and 6.8. 
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Figure 6.5 : Stress distribution of outer fiber of plate 
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Figure 6.8 : Deformed shape of plate under uniform transverse load 
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CHAPTER VII 


CONCLUSIONS 


7.1 Introduction 

The topics discussed in this chapter include the advantages of using symbolic and 
algebraic manipulation, the difficulties existing in running symbolic and algebraic software, the 
SAM in education, contributions, and the prospect of future development and application. 
Simple examples will be presented to illustrate the points where necessary. 

7.2 Advantages of symbolic and algebraic manipulation 

There are many advantages of application of symbolic and algebraic manipulation. They 
can be classified as follows : 

1. Tireless power— Together with human intelligence, the tireless capability in 
manipulating symbols and numbers has made SAM an indispensable tool in modern 
computational community. It has created the potential to challenge both previously 
intractable problems and new sophisticated formulae. Due to this advantage, the 
analytical work is pushed forward. 

2. Accuracy — A solution obtained by symbolic and algebraic manipulation is always 
exact. There is no round-off error accumulation. 

3. Reliability The resultant expressions obtained by symbolic and algebraic manipulation 
will be correct if the input information is right. In addition, the capability of automatic 
code generation eliminates any typographic errors and substantially reduces the time in 
debugging the programs. 

4. Efficiency — This is a new advantage found in this research. The symbolic template in 
nonlinear numerical analysis can significantly improve the elficiency of program 
execution. This advantage is believed to be a crucial solution in the future for fields in 
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which the development time plays an important role. For example, the real time control 
will be one of them. 

7.3 Internal swelling and mathematical limitations 

As mentioned in the last chapters, there are some difficulties existing in the symbolic 
and algebraic manipulation. The followings are the detail discussions : 

1. Memory capacity limitation — A large amount of memory space is required for symbolic 
and algebraic manipulation. This is one of its fundamental limitations. The amount of 
memory space needed for running a symbolic and algebraic manipulator varies a lot from 
one system (hardware and software) to another. Different software systems need 
different sizes of memory space, and different hardware systems may provide different 
amounts of memory space for the same software. Even the same software running in the 
same hardware system sometimes may need different memory spaces depending on 
whether external packages are connected or not. For example, 834560 bytes RAM 
(about 815 K bytes) are currently provided (Fall, 1989) for running REDUCE in the 
Michigan Terminal System (MTS) when it is invoked. If integration performance is 
involved in the computation, the external integration package should be manually 
included and the memory space will be extended to 1048560 bytes (about 1 mega bytes). 
The total memory space that MTS can provide during the computation is up to seven 
mega bytes. On the other hand, three mega bytes are provided to run REDUCE in an 
Apollo Domain workstation at the Computer Aided Engineering Network (CAEN) of 
The University of Michigan. Unlike MTS, this space can be automatically extended up 
to six mega bytes during execution if necessary. When the space requirement is beyond 
the provisions of hardware, execution will be aborted automatically. Therefore it is 
recommended that the symbolic and algebraic manipulator be implemented on a machine 
with at least one mega bytes RAM capacity to guarantee a safe execution. 

To demonstrate the mechanism of internal swell in symbolic and algebraic manipulation, 
an example is given to calculate the factorial of a number. Mathematically, a factorial is 
defined as : 

, f 1 if n = 0 ; 

\n(n - 1)! otherwise. (7.1) 

The LISP function for this problem is as follows : 
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(defun factorial (n) 


(if (= n 0) 

1 

(* n (factorial (1- n))))) 

When the function is called to calculate the factorial of four, the building process will 
occur first and then the collapsing process follows 9 . 


(factorial 4) -> 

(* 4 (factorial 3)) 

-> 

(* 4 (* 3 (factorial 2))) 

-> 

(* 4 (* 3 (* 2 (factorial 1)))) 

-> 

(* 4 (* 3 (* 2 (* 1 (factorial 0))))) 

-> 

(* 4(* 3 (* 2(* 1 1)))) 

-> 

(* 4 (* 3 (* 2 1))) 

-> 

(* 4 (* 3 2)) 

-> 

(* 4 6) 

-> 

24 


The internal swelling phenomenon occurs during the process of building. It is 
unquestionable that this phenomenon will become more serious if a larger number is 
given. Moreover if input number is negative, the recursive process will theoretically 
continue infinitely. Of course, the execution will be aborted when the provided space is 
used up. 

2. Mathematical limitation — Strictly speaking, a mathematical limitation should not be 
completely classified as limitation of symbolic and algebraic manipulation. For example, 
the analytical solution for the general 5th polynomial is proven to be non-existent. 

9 In some cases, the collapsing process may be impossible and the swelling phenomenon 
will last to the end of execution if it is not beyond the capacity of the hardware system. 
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Therefore it is also impossible for symbolic and algebraic manipulation to solve this 
problem. However, in addition to the above mentioned example, difficulties in solving 
mathematical equations whose analytical solutions exist are sometimes encountered. The 
occurrence of this phenomenon is quite system dependent. In general, most of such 
occurrences are encountered during integration and equation solving (algebraic and 
differential equation). 

7.4 Symbolic and algebraic manipulation and education 

The impact of SAM to science and engineering is significant. There are many 
publications of its application on celestial mechanics, relativity theory, and fluid mechanics. 
Compared to such successful applications, the response from educational community is far 
behind. So far, schools which officially include symbolic and algebraic manipulation in the 
content of courses include Cornell University, The University of Pennsylvania, and The 
University of Michigan. At Cornell University, MACSYMA was the first system introduced 
into the graduate-level course, in the winter of 1983. It was not until the fall of 1984 that 
Professor Richard Rand introduced the muMATH system into the sophomore engineering 
mathematics course. Unlike the MACSYMA system which ran on the mainframe, the 
muMATH system was implemented on IBM-XT and AT. At The University of Pennsylvania, 
Professor H. H. Bau employed MACSYMA in the instruction of approximate analyses. At The 
University of Michigan, Professor Noboru Kikuchi has used REDUCE to facilitate courses of 
finite element methods and applied mathematics since 1985. The other system, 
MATHEMATICA, was also implemented into the Macintosh II in the computational laboratory 
by Professor Kikuchi around 1988. Others such as Professor P. Papalambros and Professor 
R. Scott also used REDUCE in the courses on optimal design and finite element method. 

The introduction of symbolic and algebraic manipulation into the education field should 
certainly be encouraged. So far, some critic views have been reported by students at The 
University of Michigan. They are : 

1. Since there is no introductory course in symbolic and algebraic manipulation, students 
always struggle in learning the symbolic and algebraic manipulator itself rather than its 
application to the subject 

2. Most of manuals of symbolic and algebraic manipulators, such as REDUCE and 
MACSYMA, are unfriendly to the users. It is difficult for new users to understand the 
new terminologies in such a short time. 
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3. There is no appropriate textbook to facilitate to teach symbolic and algebraic 
manipulation and its application 10 . Unlike numerical analysis, the amounts of symbolic 
and algebraic manipulation results are not predictable. Therefore if assignment is not 
carefully designed, it could turn out just as simple as a symbol (say '0') or several 
hundred pages of outputs or nothing at all due to the internal swelling problem. 

4. Qualified instructors are not easily found. 

5. The software may not be fully operational. For instance, the MACSYMA system at The 
University of Michigan has just a half of its full capabilities. It is not easy to use because 
some functions cannot be found even when they are listed in the manual. The 
MATHEMATICA system is only implemented in some specific offices and is not yet 
available for public use. 

In order to overcome these problems, some proposals are suggested as follows : 

1. The education of symbolic and algebraic manipulation should start from the early 
undergraduate period. It is recommended that the existing "Numerical analysis" course 
be revised into "Numerical analysis and symbolic manipulation". The concept of 
symbolic manipulation, the use of available software, and the complementarity between 
symbolic manipulation and numerical analysis should be taught in the course. 

2. It is urgent to design a textbook for such a new course. The existing manuals need be 
revised for easy accessibility. 

3. The software systems should be rechecked and made available to the public. 

7.5 Contributions of this study 

The study presented in this report is believed to have made three original contributions 
to applied mechanics and symbolic manipulation. They are : 

A). To applied mechanics : 

1. Before this study, all of the advantages from the applications of symbolic and algebraic 
manipulation were either in handling lengthy formulae, or in increasing the accuracy of 
solution. In addition, this report points out for the first time a new advantage in 

10 The one written by Gerhand Rayna is rather an experimental book of REDUCE than an 
application of REDUCE in applied mechanics. 
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improving the efficiency of the execution of a numerical program. It is believed that this 
advantage will be crucial in such applications that the development time plays an 
important role. For instance, if the template is prepared in symbolic form beforehand and 
implemented in a chip, the data received by a heat-seeking missile can simply be 
substituted into the template. The response time will be substantially shortened. 

2. The closed-form solution of a stiffness matrix of a 4-node quadrilateral isoparametric 
element was not available before. This dissertation presents the first analytical solutions 
of it. The contributions to the finite element analysis by this breakthrough are multifolds. 
First, the integration error is eliminated and the solutions are more accurate. Secondly, 
the closed-form solution can be automatically coded into a fortran subroutine. This 
allows the element stiffness matrix to be obtained by simple substitution of nodal 
coordinates. Of course, the fortran programming is simplified and the assemblage of the 
global stiffness matrix is expedited. 

B). To Symbolic and algebraic manipulation : 

3. Although the difficulties of the internal swelling problem and mathematical limitation 
were well known in the SAM field, no one has given the remedy for it. This report 
proposes a sysmatic pre-treatment method to avoid these difficulties and then 
successfully applies it to solve the problems. 

7.6 Prospects and continuations of this research 

As the criterion of the quality of results (in both industry and academia) becomes more 

and more strict, it is expected that more and more sophisticated formulations will be produced. 

Some expectations of future trends are as follows : 

1. The design trend of symbolic and algebraic systems will continuously go towards 
smaller, more convenient packages for personal computers or even calculators. 
However, the mainframe SAM system will still co-exist to process large-expressions. 

2. Applications of SAM in industry are scare at this time. However, this situation will 
change gradually after the teaching of symbolic and algebraic manipulation is actually 
implemented in schools. 
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3. The relationship between numerical analysis and symbolic manipulation will be 
smoother in the future. The switch from symbolic manipulation to numerical analysis (or 
vice versa) is expected to be automatic eventually 

4. The gap between theoretical analysis and computational experiments will be smaller and 
smaller. 

5. The reconsideration of every problem, equation, and formulation will become necessary. 
Regardless of whether they have already been solved or not. The solved problem can be 
used to check the correctness of solutions by SAM. The unsolved problem might then 
become solvable with the employment of SAM. 

6. The inclusion of higher order terms for applied mechanics problems will become popular 
due to the availability of SAM systems. 

7. To debug the symbolic program and to check the correctness of results are the important 
associated tasks of symbolic and algebraic manipulation. It is expected that the self 
debugging function of symbolic manipulators will be developed soon. A technique for 
the sysmatic checking of results from symbolic and algebraic manipulation should be 
available in the future. 

The following three topics are closely relative to this study, and should be continued in 
future research. They are : 

1. Extension of methodology in constructing a stiffness matrix for 2-D isoparametrical 
quadrilateral element to a 3-D problem. 

2 

jL_(L) 

2. Inclusion of the higher derivative term of au ' * in Equation (6.23) to Equation (6.25). 
This was neglected in the original formulation by Lee and Kobayashi and in this thesis. 

3. Extension of methodology presented in section 7.3 to solve the general shell problem. 
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APPENDIX A 


Proof of Equation (5.12) 
Given : f { V r ) - - m * g * ^ 


Prove : f(V r ) • V] -f (P'j ■ V' r ±Q 


Proof : Let m*g =1 for simplicity, 

r r t 

( TnT' ' l/;r) 

v,v : v;‘ 

"¥J + W\ 


There are three cases for discussions : 


Case 1 : when V' r > 0 .equation (A.l) will be 

r-VUl-^Lr) [ >0 - whenV - s0 

\V r | ’ [ - 0, when V r > 0 

Case 2 : when V ' - 0 , equality is hold. 

Case 3 : when V * < 0 , equation (A.l) will be 

T - V\*( - l - Zl,) f >0 ' whenV.iO 
\VJ’ 1-0, when V r <0 

Therefore, f(V r ) • V' r -f{V' r ) ■ ?' r *0 is proven 


(A. 1) 


(A. 2) 


(A.3) 
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APPENDIX B 


Proof of Equation (5.16) 

Given : f~ ( V , ) ■ - m * g * p-j t 

Prove : jJ'V ' <&. - dV ' '< V ' ' V '~> 

Proof : Let m*g =1 for simplicity, 

rKL A < I v _ f K| 

J 0 I'd? “J 0 ~^lf Vrt 1 ~J 0 ~ sl 8 n (Vr)dV r 

ifv.>om.^r-dV'.-\v;\ ^r-dv, --|v,i 

v v , <o - jf !r • -|v;i /] p |r ■ - |v,i 

There are four cases for discussions : 

1. V ’ > 0,V r > 0 case : 


Left side of (B.l)=-|y r ’|+ |v r | — V‘ + V 
V , 

” ~ r - ^r)=rightsideof(B.l) 


2 . V r *>0,V r <0 


case : 


Left side of (B.l)=|v r *|- |V r | - + V r 

V r 

I^r( y r - y r ) =right side of (B. 1) 


s - 


3 v; <o,v r >o 


case : 


Left side of (B.l)= _ |v;| + |v r |-v; + v 


(5.1) 


(5.2) 

(5.3) 

(5.4) 
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s - v; + v , 


Vr 

F7 


( V' ’ - V r ) =right side of (B. 1) 


4 V r <0,V r <0 ca ^g ; 

Leftside of (B.l)=|v r *| - |V r | - - V’ r + V 


-(v;-v r ) = -(-T -p)(v;-v r ) 


=-l*(right side of (B. 1)) (B.5) 

Equation (B.5) implies that equality is hold and both sides of (B. 1) are zero for 

case 4. 

Therefore equation (B. 1) is proven. 
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